{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EJq4ZsOoX1qG"
      },
      "source": [
        "## LQR Example\n",
        "\n",
        "This notebook demonstrates how to solve a finite-horizon Linear Quadratic Regulator (LQR) problem using a Gaussian Factor Graph in GTSAM. The goal is to minimize a quadratic cost w.r.t. dynamics constraints over a linear dynamical system using factor-graph-based optimization. The results are compared against results from classical Riccati equations. Adapted from https://gtsam.org/2019/11/07/lqr-control.html\n",
        "\n",
        "Author(s): [Zhouyu Zhang](https://zhangzdd.github.io/zzy_webpage/)\n",
        "\n",
        "<a href=\"https://colab.research.google.com/github/borglab/gtsam/blob/develop/python/gtsam/examples/LQRExample.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bc7zas_kX3U-"
      },
      "source": [
        "GTSAM Copyright 2010-2025, Georgia Tech Research Corporation, Atlanta, Georgia 30332-0415 All Rights Reserved\n",
        "\n",
        "Authors: Frank Dellaert, et al. (see THANKS for the full author list)\n",
        "\n",
        "See LICENSE for the license information"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-wKB9H9ZaASf",
        "outputId": "7c9d3d23-34bb-4a57-9a29-f797d041e327"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Proceed (Y/n)? y\n",
            "Proceed (Y/n)? y\n"
          ]
        }
      ],
      "source": [
        "# In colab install gtsam and gtbook\n",
        "try:\n",
        "    import google.colab\n",
        "    %pip install --quiet gtsam-develop\n",
        "    %pip install --quiet gtbook\n",
        "except ImportError:\n",
        "    pass # Not in Colab\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "hUK913WGZvB1"
      },
      "outputs": [],
      "source": [
        "# Import required libraries\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import gtsam\n",
        "from gtsam import symbol\n",
        "from gtsam.symbol_shorthand import X, U\n",
        "import graphviz\n",
        "from gtbook.display import show"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zYeS1drKX938"
      },
      "source": [
        "# Discrete-Time State-Space Form (2D State: Position and Velocity)\n",
        "We look to work on a 1D double integrator system, in this sense, we define the state and input as:\n",
        "\n",
        "$$\n",
        "x_k = \\begin{bmatrix} x \\\\ \\dot{x} \\end{bmatrix}, \\quad u_k = \\ddot{x}\n",
        "$$\n",
        "\n",
        "We can think of it as a linear vehicle for which we can control the gas paddle (acceleration). Then the discrete-time dynamics become:\n",
        "\n",
        "$$\n",
        "x_{k+1} = A x_k + B u_k\n",
        "$$\n",
        "\n",
        "Where:\n",
        "\n",
        "$$\n",
        "A = \\begin{bmatrix} 1 & \\Delta t \\\\ 0 & 1 \\end{bmatrix}, \\quad\n",
        "B = \\begin{bmatrix} \\frac{1}{2} \\Delta t^2 \\\\ \\Delta t \\end{bmatrix}\n",
        "$$\n",
        "\n",
        "So:\n",
        "\n",
        "$$\n",
        "\\begin{bmatrix} x_{k+1} \\\\ \\dot{x}_{k+1} \\end{bmatrix}\n",
        "=\n",
        "\\begin{bmatrix} 1 & \\Delta t \\\\ 0 & 1 \\end{bmatrix}\n",
        "\\begin{bmatrix} x_k \\\\ \\dot{x}_k \\end{bmatrix}\n",
        "+\n",
        "\\begin{bmatrix} \\frac{1}{2} \\Delta t^2 \\\\ \\Delta t \\end{bmatrix} u_k\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AM-cyKjEYBut"
      },
      "source": [
        "# LQR Problem Setup\n",
        "\n",
        "We seek to minimize the infinite-horizon cost function:\n",
        "\n",
        "$$\n",
        "J = \\sum_{k=0}^\\infty \\left( x_k^\\top Q x_k + u_k^\\top R u_k \\right)\n",
        "$$\n",
        "\n",
        "Choose the weighting matrices:\n",
        "\n",
        "$$\n",
        "Q = \\begin{bmatrix} q_1 & 0 \\\\ 0 & q_2 \\end{bmatrix}, \\quad\n",
        "R = r\n",
        "$$\n",
        "\n",
        "- $ Q \\succeq 0 $: penalizes position and velocity\n",
        "- $ R \\succ 0 $: penalizes control effort\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "iQRHn0qFdYXL"
      },
      "outputs": [],
      "source": [
        "# Define system dynamics: x_{k+1} = A x_k + B u_k\n",
        "A = np.array([[1.0, 0.1],\n",
        "        [0.0, 1.0]])\n",
        "B = np.array([[0.005],\n",
        "        [1.0]])\n",
        "\n",
        "# Define cost matrices: minimize sum of (xᵀQx + uᵀRu)\n",
        "Q = np.eye(2)\n",
        "R = np.array([[1.0]])\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WvdqMTD5ZzmV"
      },
      "source": [
        "\n",
        "###  Graphical Model Representation of Finite-Horizon LQR\n",
        "\n",
        "We formulate the finite-horizon Linear Quadratic Regulator (LQR) problem as a **factor graph**, allowing us to interpret optimal control as a **maximum a posteriori (MAP)** estimation problem.\n",
        "\n",
        "---\n",
        "\n",
        "###  Problem Statement\n",
        "\n",
        "Given discrete-time linear dynamics:\n",
        "\n",
        "$$\n",
        "x_{k+1} = A x_k + B u_k\n",
        "$$\n",
        "\n",
        "and a quadratic cost function:\n",
        "\n",
        "$$\n",
        "J(x, u) = \\sum_{k=0}^{N-1} \\left( x_k^\\top Q x_k + u_k^\\top R u_k \\right) + x_N^\\top Q x_N\n",
        "$$\n",
        "\n",
        "we aim to find the state-control sequence $ \\{x_k, u_k\\}_{k=0}^{N-1} $ minimizing $ J $, subject to initial state $ x_0 $ and dynamics.\n",
        "\n",
        "---\n",
        "\n",
        "### MAP Estimation View\n",
        "\n",
        "This control problem is equivalent to a **MAP estimation** over a factor graph:\n",
        "\n",
        "$$\n",
        "\\arg\\max_{x_{0:N}, u_{0:N-1}} \\; p(x_{0:N}, u_{0:N-1} \\mid \\text{costs, dynamics})\n",
        "$$\n",
        "\n",
        "This posterior can be written as a product of factors:\n",
        "\n",
        "$$\n",
        "p(x, u) \\propto \\phi_0(x_0)\n",
        "\\prod_{k=0}^{N-1} \\phi_Q(x_k) \\phi_R(u_k) \\phi_f(x_k, u_k, x_{k+1})\n",
        "\\cdot \\phi_Q(x_N)\n",
        "$$\n",
        "\n",
        "---\n",
        "\n",
        "###  Prior Factor (Initial Condition)\n",
        "\n",
        "We encode the prior on the initial state as a Gaussian:\n",
        "\n",
        "$$\n",
        "\\phi_0(x_0) \\propto \\exp\\left( -\\frac{1}{2} \\| x_0 - x_{\\text{init}} \\|^2_{\\Sigma_0} \\right)\n",
        "$$\n",
        "\n",
        "where $ \\Sigma_0 \\approx 0 \\cdot I $ enforces the initial condition tightly.\n",
        "\n",
        "---\n",
        "\n",
        "###  Quadratic Cost Factors\n",
        "\n",
        "Each quadratic cost $ x_k^\\top Q x_k $ is rewritten via Cholesky decomposition $ Q = L_x^\\top L_x $ as:\n",
        "\n",
        "$$\n",
        "x_k^\\top Q x_k = \\| L_x x_k \\|^2 \\quad \\Rightarrow \\quad\n",
        "\\phi_Q(x_k) \\propto \\exp\\left( -\\frac{1}{2} \\| L_x x_k \\|^2 \\right)\n",
        "$$\n",
        "\n",
        "Similarly, for control:\n",
        "\n",
        "$$\n",
        "\\phi_R(u_k) \\propto \\exp\\left( -\\frac{1}{2} \\| L_u u_k \\|^2 \\right), \\quad R = L_u^\\top L_u\n",
        "$$\n",
        "\n",
        "These are **soft constraints**, penalizing deviation from zero.\n",
        "\n",
        "---\n",
        "\n",
        "###  Dynamics Constraints (Hard Constraints)\n",
        "\n",
        "Each dynamics equation:\n",
        "\n",
        "$$\n",
        "x_{k+1} = A x_k + B u_k\n",
        "$$\n",
        "\n",
        "is encoded as a **hard constraint**, modeled as a Gaussian with **zero variance**:\n",
        "\n",
        "$$\n",
        "\\phi_f(x_k, u_k, x_{k+1}) \\propto \\delta(x_{k+1} - A x_k - B u_k)\n",
        "= \\lim_{\\Sigma \\to 0} \\exp\\left( -\\frac{1}{2} \\| x_{k+1} - A x_k - B u_k \\|^2_{\\Sigma} \\right)\n",
        "$$\n",
        "\n",
        "This enforces the equality exactly during optimization.\n",
        "\n",
        "---\n",
        "\n",
        "###  Terminal Cost\n",
        "\n",
        "The terminal cost is also encoded as a soft quadratic factor:\n",
        "\n",
        "$$\n",
        "\\phi_Q(x_N) \\propto \\exp\\left( -\\frac{1}{2} \\| L_x x_N \\|^2 \\right)\n",
        "$$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "R6YXu9dcTkWG"
      },
      "outputs": [],
      "source": [
        "# Horizon of the problem (i.e. how many discrete timesteps we want to consider)\n",
        "# Note that in python the initial index starts at 0, so the terminal states would be N-1\n",
        "N = 4\n",
        "# Cost factors: L such that Q = LᵀL\n",
        "Lx = np.linalg.cholesky(Q)\n",
        "Lu = np.linalg.cholesky(R)\n",
        "\n",
        "# Create a Gaussian factor graph for the linear-quadratic problem\n",
        "graph = gtsam.GaussianFactorGraph()\n",
        "initial_values = gtsam.Values()\n",
        "\n",
        "# Set initial state: x_0 = [0, 1]^T, add a strong prior factor to fix the initial state\n",
        "x0 = np.array([0.0, 1.0])\n",
        "prior_noise = gtsam.noiseModel.Isotropic.Sigma(2, 1e-20)\n",
        "graph.add(gtsam.JacobianFactor(X(0), np.eye(2), x0, prior_noise))\n",
        "\n",
        "# Cost + dynamics\n",
        "for k in range(N):\n",
        "    xk = X(k)\n",
        "    uk = U(k)\n",
        "\n",
        "    # Cost factor: x_k^T Q x_k\n",
        "    graph.add(gtsam.JacobianFactor(xk, Lx, np.zeros(2), gtsam.noiseModel.Unit.Create(2)))\n",
        "\n",
        "    # Cost factor: u_k^T R u_k, note that there is no cost term for terminal input\n",
        "    if k < N-1:\n",
        "        graph.add(gtsam.JacobianFactor(uk, Lu, np.zeros(1), gtsam.noiseModel.Unit.Create(1)))\n",
        "\n",
        "    # Dynamics hard constraint factor: x_{k+1} = A x_k + B u_k\n",
        "    if k <= N-2:\n",
        "        xk1 = X(k+1)\n",
        "        graph.add(gtsam.JacobianFactor(\n",
        "            xk1, np.eye(2),\n",
        "            xk, -A,\n",
        "            uk, -B,\n",
        "            np.zeros((2, 1)),\n",
        "            gtsam.noiseModel.Diagonal.Sigmas(np.zeros(2))  # hard constraint\n",
        "        ))\n",
        "\n",
        "# Terminal cost: x_N-1^T Q x_N-1\n",
        "xN = X(N-1)\n",
        "graph.add(gtsam.JacobianFactor(xN, Lx, np.zeros(2), gtsam.noiseModel.Unit.Create(2)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 464
        },
        "id": "Pg5n_cfVd5V7",
        "outputId": "1804024e-e11a-4191-d296-050a9a7fd3cb"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[0.        , 1.        ],\n",
              "       [0.0968027 , 0.36053946],\n",
              "       [0.13167395, 0.12400146],\n",
              "       [0.14365837, 0.04085496]])"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAF2CAYAAACmtO2KAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYVtJREFUeJzt3Xd8FHX+x/HX7ibZJIQklJAEiIQmvXg0AenxEAHFRvOkiCKIhR96p5wniHeKnKdyniiICiqiFBUVVIRQBERBEEEFpIQikNATSM/u/P7YZJMlhQSSTMr7+XjMg+zMd2Y/m9mQfef7ne9YDMMwEBERERERqUCsZhcgIiIiIiJS3BR0RERERESkwlHQERERERGRCkdBR0REREREKhwFHRERERERqXAUdEREREREpMJR0BERERERkQpHQUdERERERCocBR0REREREalwFHRERCqJUaNGERkZacpzWywWnnnmGVOeu6JZt24dFouFdevWmV3KFenZsyc9e/Y0uwwRqQQUdESkzIuJieGhhx7i2muvxd/fH39/f5o3b86ECRPYuXMn4PrwZLFYLrtkfdhOS0vjv//9L9dddx2BgYEEBwfTokULxo4dy549e/KsY/fu3VgsFnx9fTl//nyebbLqaNy4cZ7bV61a5a5l6dKlV/29qWhOnTrFo48+StOmTfHz86NWrVp07NiRJ554gosXL7rbLVy4kJkzZ17x8yQlJfHMM88Ua1gYNWpUod6Do0aNKrbnLAm//fYbzzzzDIcOHTK7FBGRq+JldgEiIgVZvnw5Q4YMwcvLi7vvvps2bdpgtVrZs2cPn3zyCW+88QYxMTE89dRT3Hfffe79tm7dyquvvsrf//53mjVr5l7funVrAO644w6++uorhg0bxv333096ejp79uxh+fLldOnShaZNm+aqZcGCBYSFhXHu3DmWLl3q8Xw5+fr6sn//frZs2ULHjh09tn3wwQf4+vqSkpJSHN+eCuXs2bO0b9+ehIQE7r33Xpo2bcqZM2fYuXMnb7zxBuPHjycgIABwBZ1ffvmFiRMnXtFzJSUlMW3aNIBi61144IEHiIqKcj+OiYlhypQpjB07lm7durnXN2zY8Kqep3v37iQnJ+Pj43NVx8nPb7/9xrRp0+jZs2eJ9AB+8803xX5MEZG8KOiISJl14MABhg4dSr169YiOjiY8PNxj+4wZM3j99dexWq3ceOONHtt8fX159dVXufHGG3N9kN26dSvLly/nueee4+9//7vHttdeey3P3hrDMFi4cCHDhw8nJiaGDz74IN+g07BhQzIyMvjwww89gk5KSgqffvop/fv35+OPPy7Cd6JyePvttzly5AibNm2iS5cuHtsSEhJK7IN9cencuTOdO3d2P/7xxx+ZMmUKnTt35i9/+Uu++yUmJlKlSpVCP4/VasXX1/eqajVDUlIS/v7+Zf48ikjFoaFrIlJm/fvf/yYxMZF58+blCjkAXl5ePPLII0RERBTpuAcOHACga9euubbZbDZq1KiRa/2mTZs4dOgQQ4cOZejQoXz77bf88ccf+T7HsGHDWLRoEU6n073uiy++ICkpicGDBxe61tjYWEaPHk3dunWx2+2Eh4dz66235hpW9Prrr9OiRQvsdju1a9dmwoQJ+Q6vA0hPT6d69eqMHj0617aEhAR8fX15/PHH3etSU1OZOnUqjRo1wm63ExERwd/+9jdSU1M99k1NTeX//u//CAkJoWrVqtxyyy0Ffp9yOnDgADabjeuvvz7XtsDAQPeH+549e7JixQoOHz7sHg6W1fOQlpbGlClTaNeuHUFBQVSpUoVu3bqxdu1a97EOHTpESEgIANOmTcs1rBFgz5493HnnnVSvXh1fX1/at2/P559/XqjXUZD58+djsVhYv349Dz74ILVq1aJu3boAHD58mAcffJAmTZrg5+dHjRo1uOuuu3Kd6/yu0fnhhx+46aabCAoKwt/fnx49erBp06ZcNRw7dowxY8ZQu3Zt7HY79evXZ/z48aSlpTF//nzuuusuAHr16uX+3uR8rsK813r27EnLli3Ztm0b3bt3x9/f3/1Hhbyu0Sns+2vVqlXccMMNBAcHExAQQJMmTXL9sUJEJIt6dESkzFq+fDmNGjWiU6dOxXrcevXqAa5hZF27dsXL6/L/FX7wwQc0bNiQDh060LJlS/z9/fnwww/561//mmf74cOHu68B6d27N+AabtWnTx9q1apV6FrvuOMOfv31Vx5++GEiIyM5efIkq1at4siRI+4P98888wzTpk0jKiqK8ePHs3fvXt544w22bt3Kpk2b8Pb2znVcb29vbrvtNj755BPmzJnj8Vf2ZcuWkZqaytChQwFwOp3ccsstbNy4kbFjx9KsWTN27drFK6+8wu+//86yZcvc+953330sWLCA4cOH06VLF9asWUP//v0L9Vrr1auHw+Hg/fffZ+TIkfm2e+qpp4iPj+ePP/7glVdeAXAPaUtISOCtt95yD0m8cOECb7/9Nn379mXLli20bduWkJAQ91C42267jdtvvx3IHtb466+/0rVrV+rUqcOTTz5JlSpVWLx4MYMGDeLjjz/mtttuK9TrKciDDz5ISEgIU6ZMITExEXD1NH733XcMHTqUunXrcujQId544w169uzJb7/9hr+/f77HW7NmDf369aNdu3ZMnToVq9XKvHnz6N27Nxs2bHD3LB4/fpyOHTty/vx5xo4dS9OmTTl27BhLly4lKSmJ7t2788gjj+Qa9pn1b1Hea2fOnKFfv34MHTqUv/zlL4SGhuZZe2HfX7/++isDBgygdevWPPvss9jtdvbv359nmBMRAcAQESmD4uPjDcAYNGhQrm3nzp0zTp065V6SkpJytVmyZIkBGGvXrs21zel0Gj169DAAIzQ01Bg2bJgxa9Ys4/Dhw3nWkpaWZtSoUcN46qmn3OuGDx9utGnTJlfbHj16GC1atDAMwzDat29vjBkzxl2zj4+P8e677xpr1641AGPJkiUFfg/OnTtnAMaLL76Yb5uTJ08aPj4+xp///GfD4XC417/22msGYLzzzjvudSNHjjTq1avnfrxy5UoDML744guPY958881GgwYN3I/ff/99w2q1Ghs2bPBoN3v2bAMwNm3aZBiGYezYscMAjAcffNCj3fDhww3AmDp1aoGvNzY21ggJCTEAo2nTpsa4ceOMhQsXGufPn8/Vtn///h6vJUtGRoaRmprqse7cuXNGaGioce+997rXnTp1Kt+a+vTpY7Rq1cpISUlxr3M6nUaXLl2Mxo0bF/gactq6dasBGPPmzXOvmzdvngEYN9xwg5GRkeHRPq/38ebNmw3AeO+999zrst4/We9tp9NpNG7c2Ojbt6/hdDo9jle/fn3jxhtvdK8bMWKEYbVaja1bt+Z6rqx98/vZKcp7Levna/bs2bmep0ePHkaPHj3cjwv7/nrllVcMwDh16lSuY4qI5EVD10SkTEpISACy/1KfU8+ePQkJCXEvs2bNKtKxLRYLK1eu5F//+hfVqlXjww8/ZMKECdSrV48hQ4bkGobz1VdfcebMGYYNG+ZeN2zYMH7++Wd+/fXXfJ9n+PDhfPLJJ6SlpbF06VJsNluRegP8/Pzw8fFh3bp1nDt3Ls82q1evJi0tjYkTJ2K1Zv+Xfv/99xMYGMiKFSvyPX7v3r2pWbMmixYtcq87d+4cq1atYsiQIe51S5YsoVmzZjRt2pTTp0+7l6yeqqxhYV9++SUAjzzyiMfzFHbCgNDQUH7++WfGjRvHuXPnmD17NsOHD6dWrVr885//xDCMyx7DZrO5e6ecTidnz54lIyOD9u3bs3379svuf/bsWdasWcPgwYO5cOGC+7WeOXOGvn37sm/fPo4dO1ao11OQ+++/H5vN5rHOz8/P/XV6ejpnzpyhUaNGBAcHF1j7jh072LdvH8OHD+fMmTPumhMTE+nTpw/ffvstTqcTp9PJsmXLGDhwIO3bt891HIvFUmDNRX2v2e32PIdGXqqw76/g4GAAPvvsM48hoSIi+VHQEZEyqWrVqgAeUwpnmTNnDqtWrWLBggVXfHy73c5TTz3F7t27OX78OB9++CHXX389ixcv5qGHHvJou2DBAurXr+8eKrN//34aNmyIv78/H3zwQb7PMXToUOLj4/nqq6/44IMPGDBggPt15ZSWlkZsbKzH4nA4sNvtzJgxg6+++orQ0FC6d+/Ov//9b2JjY937Hj58GIAmTZp4HNPHx4cGDRq4t+fFy8uLO+64g88++8x9LcQnn3xCenq6R9DZt28fv/76q0e4DAkJ4dprrwXg5MmT7lqsVmuuWcUura0g4eHhvPHGG5w4cYK9e/fy6quvuod4vf3224U6xrvvvkvr1q3x9fWlRo0ahISEsGLFCuLj4y+77/79+zEMg6effjrX6506darH670a9evXz7UuOTmZKVOmEBERgd1up2bNmoSEhHD+/PkCa9+3bx8AI0eOzFXzW2+9RWpqKvHx8Zw6dYqEhARatmx5RTUX9b1Wp06dQk08UNj315AhQ+jatSv33XcfoaGhDB06lMWLFyv0iEi+dI2OiJRJQUFBhIeH88svv+TalnXNTnHd5yM8PJyhQ4dyxx130KJFCxYvXsz8+fPx8vIiISGBL774gpSUlDzvjbNw4UKee+65PP8aHh4eTs+ePXnppZfYtGlTvjOtfffdd/Tq1ctjXUxMDJGRkUycOJGBAweybNkyVq5cydNPP8306dNZs2YN11133VW/9qFDhzJnzhy++uorBg0axOLFi2natClt2rRxt3E6nbRq1YqXX345z2MUdTKIwrBYLFx77bVce+219O/fn8aNGxc4012WBQsWMGrUKAYNGsRf//pXatWqhc1mY/r06e5JKAqS9aH58ccfp2/fvnm2adSoUdFf0CVy9t5kefjhh5k3bx4TJ06kc+fOBAUFYbFYGDp0aIEf5rO2vfjii7Rt2zbPNgEBAZw9e/aq6y6KvF5jXgr7/vLz8+Pbb79l7dq1rFixgq+//ppFixbRu3dvvvnmm1w9ZCIiCjoiUmb179+ft956K8/70ZQEb29vWrduzb59+zh9+jRhYWF88sknpKSk8MYbb1CzZk2P9nv37uUf//gHmzZt4oYbbsjzmMOHD+e+++4jODiYm2++Oc82bdq0YdWqVR7rwsLC3F83bNiQxx57jMcee4x9+/bRtm1bXnrpJRYsWOCeWGHv3r00aNDAvU9aWhoxMTEe93XJS/fu3QkPD2fRokXccMMNrFmzhqeeesqjTcOGDfn555/p06dPgcOb6tWrh9Pp5MCBAx5/9d+7d2+BNVxOgwYNqFatGidOnHCvy6+OpUuX0qBBAz755BOPNlm9MZfbP+t76O3tfdnvXXFbunQpI0eO5KWXXnKvS0lJKXD2PMi+L09gYGCBNYeEhBAYGJjnHw9yyu97c7XvtfwU9v0Frqm1+/TpQ58+fXj55Zd5/vnneeqpp1i7dm2pny8RKfs0dE1Eyqy//e1v+Pv7c++99xIXF5dre2Gu2cjLvn37OHLkSK7158+fZ/PmzVSrVs09/fCCBQto0KAB48aN48477/RYHn/8cQICAgocvnbnnXcydepUXn/99XyH8VSrVo2oqCiPxdfXl6SkpFw3Fm3YsCFVq1Z1DzWLiorCx8eHV1991eP78fbbbxMfH3/ZGc+sVit33nknX3zxBe+//z4ZGRkew9YABg8ezLFjx5g7d26u/ZOTk92zhvXr1w+AV1991aPNzJkzC6whyw8//OA+Vk5btmzhzJkzHuGpSpUqeQ7nyvqrfs7vxQ8//MDmzZs92mXNYHZpiKhVqxY9e/Zkzpw5HsEqy6lTpwr1Wq6EzWbL9Z7+3//+h8PhKHC/du3a0bBhQ/7zn//kOdQzq2ar1cqgQYP44osv+PHHH3O1y3rurHv6XPq9udr3Wn4K+/7Kq0cqqwfr0mmoRURAPToiUoY1btyYhQsXMmzYMJo0acLdd99NmzZtMAyDmJgYFi5ciNVqdd+HpLB+/vlnhg8fTr9+/ejWrRvVq1fn2LFjvPvuuxw/fpyZM2dis9k4fvw4a9euzXVxfRa73U7fvn1ZsmQJr776ap7TOAcFBXncn6Uofv/9d/r06cPgwYNp3rw5Xl5efPrpp8TFxbmnfg4JCWHy5MlMmzaNm266iVtuuYW9e/fy+uuv06FDhwJvVJllyJAh/O9//2Pq1Km0atXKPZVwlnvuuYfFixczbtw41q5dS9euXXE4HOzZs4fFixezcuVK2rdvT9u2bRk2bBivv/468fHxdOnShejoaPbv31+o1/v+++/zwQcfcNttt9GuXTt8fHzYvXs377zzDr6+vh73S2nXrh2LFi1i0qRJdOjQgYCAAAYOHMiAAQP45JNPuO222+jfvz8xMTHMnj2b5s2be4QAPz8/mjdvzqJFi7j22mupXr06LVu2pGXLlsyaNYsbbriBVq1acf/999OgQQPi4uLYvHkzf/zxBz///HOhXk9RDRgwgPfff5+goCCaN2/O5s2bWb16dZ73dcrJarXy1ltv0a9fP1q0aMHo0aOpU6cOx44dY+3atQQGBvLFF18A8Pzzz/PNN9/Qo0cP91TOJ06cYMmSJWzcuJHg4GDatm2LzWZjxowZxMfHY7fb6d27N7Vq1brq91peCvv+evbZZ/n222/p378/9erV4+TJk7z++uvUrVs33x5VEankTJvvTUSkkPbv32+MHz/eaNSokeHr62v4+fm5px/esWNHnvsUNL10XFyc8cILLxg9evQwwsPDDS8vL6NatWpG7969jaVLl7rbvfTSSwZgREdH51vb/PnzDcD47LPPDMPwnF46P4WdXvr06dPGhAkTjKZNmxpVqlQxgoKCjE6dOhmLFy/O1fa1114zmjZtanh7exuhoaHG+PHjjXPnznm0uXR66SxOp9OIiIgwAONf//pXnrWkpaUZM2bMMFq0aGHY7XajWrVqRrt27Yxp06YZ8fHx7nbJycnGI488YtSoUcOoUqWKMXDgQOPo0aOFml56586dxl//+lfjT3/6k1G9enXDy8vLCA8PN+666y5j+/btHm0vXrxoDB8+3AgODjYA9+tyOp3G888/b9SrV8+w2+3GddddZyxfvjzP1/7dd98Z7dq1M3x8fHLVd+DAAWPEiBFGWFiY4e3tbdSpU8cYMGCAx/vjcgqaXjqv6Z3PnTtnjB492qhZs6YREBBg9O3b19izZ49Rr149Y+TIke52l04vneWnn34ybr/9dqNGjRqG3W436tWrZwwePDjX+/fw4cPGiBEjjJCQEMNutxsNGjQwJkyY4DEt99y5c40GDRoYNpst13MV5r1W0M/BpdNLG0bh3l/R0dHGrbfeatSuXdvw8fExateubQwbNsz4/fff83weERGLYVzh2A8REREpddHR0URFRbFhwwb1ZIiIFEDX6IiIiJQjWdcOXTo5hoiIeFKPjoiISDmQmJjIBx98wH//+18SEhLc9y0SEZG86X9IERGRcuDUqVM8/PDD+Pn58fHHHyvkiIhchnp0RERERESkwtGfg0REREREpMJR0BERERERkQqnXNww1Ol0cvz4capWrYrFYjG7HBERERERMYlhGFy4cIHatWsXeL1iuQg6x48fJyIiwuwyRERERESkjDh69Ch169bNd3u5CDpVq1YFXC8mMDDQ5GpERERERMQsCQkJREREuDNCfspF0MkarhYYGKigIyIiIiIil72kRZMRiIiIiIhIhaOgIyIiIiIiFY6CjoiIiIiIVDjl4hodEREREZGS4HA4SE9PN7sMycHb2xubzXbVx1HQEREREZFKxzAMYmNjOX/+vNmlSB6Cg4MJCwu7qntoKuiIiIiISKWTFXJq1aqFv7+/bkpfRhiGQVJSEidPngQgPDz8io+loCMiIiIilYrD4XCHnBo1aphdjlzCz88PgJMnT1KrVq0rHsZW5MkIvv32WwYOHEjt2rWxWCwsW7bssvusW7eOP/3pT9jtdho1asT8+fOvoFQRERERkauXdU2Ov7+/yZVIfrLOzdVcP1XkoJOYmEibNm2YNWtWodrHxMTQv39/evXqxY4dO5g4cSL33XcfK1euLHKxIiIiIiLFRcPVyq7iODdFHrrWr18/+vXrV+j2s2fPpn79+rz00ksANGvWjI0bN/LKK6/Qt2/foj69+dKTwdvP7CpERERERKQAJX4fnc2bNxMVFeWxrm/fvmzevDnffVJTU0lISPBYygRHBrx/OyybAGlJZlcjIiIiIlIs1q1bh8ViuewsdJGRkcycObNUarpaJR50YmNjCQ0N9VgXGhpKQkICycnJee4zffp0goKC3EtERERJl1k4hzfB0e9hxwJ4KwpO7ze7IhERERGpREaNGoXFYsFiseDj40OjRo149tlnycjIuKrjdunShRMnThAUFATA/PnzCQ4OztVu69atjB079qqeq7SUeNC5EpMnTyY+Pt69HD161OySXBr0gHuWQZUQOPkrvNkTfv3U7KpEREREpBK56aabOHHiBPv27eOxxx7jmWee4cUXX7yqY/r4+BTqvjUhISHlZhKHEg86YWFhxMXFeayLi4sjMDDQPXXcpex2O4GBgR5LmdGgBzywAa7pAmkXYMko+OoJyEgzuzIRERERqQTsdjthYWHUq1eP8ePHExUVxeeff865c+cYMWIE1apVw9/fn379+rFv3z73focPH2bgwIFUq1aNKlWq0KJFC7788kvAc+jaunXrGD16NPHx8e7eo2eeeQbIPXTtyJEj3HrrrQQEBBAYGMjgwYM9Pvs/88wztG3blvfff5/IyEiCgoIYOnQoFy5cKPHvU4kHnc6dOxMdHe2xbtWqVXTu3Lmkn7rkBIbDyC+g60TX4x9mw7x+cL6M9DyJiIiISJEYhkFSWoYpi2EYV1W7n58faWlpjBo1ih9//JHPP/+czZs3YxgGN998s3uK5gkTJpCamsq3337Lrl27mDFjBgEBAbmO16VLF2bOnElgYCAnTpzgxIkTPP7447naOZ1Obr31Vs6ePcv69etZtWoVBw8eZMiQIR7tDhw4wLJly1i+fDnLly9n/fr1vPDCC1f1mgujyLOuXbx4kf37s69NiYmJYceOHVSvXp1rrrmGyZMnc+zYMd577z0Axo0bx2uvvcbf/vY37r33XtasWcPixYtZsWJF8b0KM9i84MZpcM318OkDcOxHmNMNbp8LjW80uzoRERERKYLkdAfNp5hz+5Pfnu2Lv0+RP5ZjGAbR0dGsXLmSfv36sWzZMjZt2kSXLl0A+OCDD4iIiGDZsmXcddddHDlyhDvuuINWrVoB0KBBgzyP6+PjQ1BQEBaLhbCwsHyfPzo6ml27dhETE+O+pv69996jRYsWbN26lQ4dOgCuQDR//nyqVq0KwD333EN0dDTPPfdckV9zURS5R+fHH3/kuuuu47rrrgNg0qRJXHfddUyZMgWAEydOcOTIEXf7+vXrs2LFClatWkWbNm146aWXeOutt8rn1NJ5adIPHvgWwttC8jn44E6IftY1Q5uIiIiISDFbvnw5AQEB+Pr60q9fP4YMGcKoUaPw8vKiU6dO7nY1atSgSZMm7N69G4BHHnmEf/3rX3Tt2pWpU6eyc+fOq6pj9+7dREREeEwc1rx5c4KDg93PCa7hblkhByA8PJyTJ09e1XMXRpGjY8+ePQvsXps/f36e+/z0009Ffaryo1okjPkGVv4dtr4FG16Co1vgjrehauhldxcRERERc/l52/jtWXP+EO/nbStS+169evHGG2/g4+ND7dq18fLy4vPPP7/sfvfddx99+/ZlxYoVfPPNN0yfPp2XXnqJhx9++EpLLxRvb2+PxxaLBafTWaLPCWV01rVyycsO/V9yhRvvKnBog2so26GNZlcmIiIiIpdhsVjw9/EyZbncTGeXqlKlCo0aNeKaa67By8vVb9GsWTMyMjL44Ycf3O3OnDnD3r17ad68uXtdREQE48aN45NPPuGxxx5j7ty5eT6Hj48PDoejwDqaNWvG0aNHPWZI/u233zh//rzHc5pFQae4tboTxq6DkGZwMQ7eHQgbXoZSSK0iIiIiUjk1btyYW2+9lfvvv5+NGzfy888/85e//IU6depw6623AjBx4kRWrlxJTEwM27dvZ+3atTRr1izP40VGRnLx4kWio6M5ffo0SUlJudpERUXRqlUr7r77brZv386WLVsYMWIEPXr0oH379iX6egtDQackhFwL90dD66FgOCF6Gnw0DJLOml2ZiIiIiFRQ8+bNo127dgwYMIDOnTtjGAZffvmle+iYw+FgwoQJNGvWjJtuuolrr72W119/Pc9jdenShXHjxjFkyBBCQkL497//nauNxWLhs88+o1q1anTv3p2oqCgaNGjAokWLSvR1FpbFuNr57EpBQkICQUFBxMfHl6176lyOYcD2d+HLv4EjFYKugcHzoU47sysTERERqbRSUlKIiYmhfv36+Pr6ml2O5KGgc1TYbKAenZJksUC7UXDfKqhWH+KPwDs3wZa5rhAkIiIiIiIlQkGnNIS3gQfWQ9MB4EiDLx+Hj8dAasnfEVZEREREpDJS0CktvkEwZAH8+TmwesEvH8ObvSDuN7MrExERERGpcBR0SpPFAl0eglFfQtXacGYfzO0NOxaaXZmIiIiISIWioGOGazrBuA3QsDdkJMOy8fD5w5CebHZlIiIiIiIVgoKOWarUhLuXQs+/AxbY/h68fSOcOWB2ZSIiIiIi5Z6CjpmsNuj5BNzzCfjXhNhd8GZP+O1zsysTERERESnXFHTKgoa9XUPZIq6H1ARYfA98PRky0syuTERERESkXFLQKSsCa8Oo5dDlYdfj71+H+f0h/g9z6xIRERERKYcUdMoSmzf8+V8w5AOwB8EfW2B2N9i/2uzKRERERKQCiIyMZObMmWX2eMVJQacsajbAdYPRsNaQfBYW3AlrngOnw+zKRERERMQkAwcO5Kabbspz24YNG7BYLOzcubNUa9q6dStjx451P7ZYLCxbtqxUa8iPgk5ZVb0+jFkF7UYDBnz7b3j/Nrh4yuzKRERERMQEY8aMYdWqVfzxR+5LG+bNm0f79u1p3bp1qdYUEhKCv79/qT5nYSnolGXevjBwJtw+F7z9IWY9zL4BDn9ndmUiIiIiUsoGDBhASEgI8+fP91h/8eJFlixZwpgxY9i4cSPdunXDz8+PiIgIHnnkERITE/M95pEjR7j11lsJCAggMDCQwYMHExcX59Hmiy++oEOHDvj6+lKzZk1uu+0297acQ9ciIyMBuO2227BYLERGRnLo0CGsVis//vijxzFnzpxJvXr1cDqdV/4NuQwFnfKg9WC4fy3UbAIXY2H+ANj0XzAMsysTERERqRgMA9ISzVkK+ZnOy8uLESNGMH/+fIwc+yxZsgSHw0Hnzp256aabuOOOO9i5cyeLFi1i48aNPPTQQ3kez+l0cuutt3L27FnWr1/PqlWrOHjwIEOGDHG3WbFiBbfddhs333wzP/30E9HR0XTs2DHP423duhVw9S6dOHGCrVu3EhkZSVRUFPPmzfNoO2/ePEaNGoXVWnJxxGIYZf/TckJCAkFBQcTHxxMYGGh2OeZJvQjLJ8KuJa7HTW6GQa+DXzVTyxIREREpT1JSUoiJiaF+/fr4+vq6VqYlwvO1zSno78fBp0qhmu7Zs4dmzZqxdu1aevbsCUD37t2pV68edrsdm83GnDlz3O03btxIjx49SExMxNfXl8jISCZOnMjEiRNZtWoV/fr1IyYmhoiICAB+++03WrRowZYtW+jQoQNdunShQYMGLFiwIM96ch4PXNfofPrppwwaNMjdZvHixYwbN44TJ05gt9vZvn077du35+DBg+5eoEvleY4yFTYbqEenPLEHuIax9X8ZbD6w90uY0wOO/2R2ZSIiIiJSCpo2bUqXLl145513ANi/fz8bNmxgzJgx/Pzzz8yfP5+AgAD30rdvX5xOJzExMbmOtXv3biIiItwhB6B58+YEBweze/duAHbs2EGfPn2uquZBgwZhs9n49NNPAZg/fz69evXKN+QUF68SPboUP4sFOoyBOn+CxSPh/GF4+89w03RoP8a1XURERESKxtvf1bNi1nMXwZgxY3j44YeZNWsW8+bNo2HDhvTo0YOLFy/ywAMP8Mgjj+Ta55prrrmi0vz8/K5ov5x8fHwYMWIE8+bN4/bbb2fhwoX897//verjXo6CTnlV+zrXFNTLHnT17Kx4DI58DwNmunp+RERERKTwLJZCDx8z2+DBg3n00UdZuHAh7733HuPHj8disfCnP/2J3377jUaNGhXqOM2aNePo0aMcPXrUY+ja+fPnad68OQCtW7cmOjqa0aNHF+qY3t7eOBy5b4ly33330bJlS15//XUyMjK4/fbbC/lqr5yGrpVnftVg6EK48Vmw2FzX7sztDSf3mF2ZiIiIiJSQgIAAhgwZwuTJkzlx4gSjRo0C4IknnuC7777joYceYseOHezbt4/PPvss38kIoqKiaNWqFXfffTfbt29ny5YtjBgxgh49etC+fXsApk6dyocffsjUqVPZvXs3u3btYsaMGfnWFhkZSXR0NLGxsZw7d869vlmzZlx//fU88cQTDBs2rFh6ii5HQae8s1ig66MwajlUDYfTe2FuL/h5kdmViYiIiEgJGTNmDOfOnaNv377Uru2aRKF169asX7+e33//nW7dunHdddcxZcoU9/ZLWSwWPvvsM6pVq0b37t2JioqiQYMGLFqU/TmyZ8+eLFmyhM8//5y2bdvSu3dvtmzZkm9dL730EqtWrSIiIoLrrrsuV81paWnce++9xfAduDzNulaRXDwFH49x3W8HoN0ouGmG6348IiIiIgIUPKOXlJx//vOfLFmyhJ07d162rWZdE08BIXDPp9DjCcAC2+bD2zfC2YNmVyYiIiIildTFixf55ZdfeO2113j44YdL7XkVdCoaqw16/R3+8jH4VYfYnTCnJ+xebnZlIiIiIlIJPfTQQ7Rr146ePXuW2rA1UNCpuBr1gXEboG5HSI2HRXfDyqfAkW52ZSIiIiJSicyfP5/U1FQWLVqEzWYrtedV0KnIgurC6C/h+gmux5tfg/kDIMGkOeJFREREREqJgk5FZ/OGm56Hwe+DPRCOfg+zb4ADa8yuTERERESkxCjoVBbNb4Gx6yCsFSSdgfdvh3UvgDP3DZ1EREREKgOn02l2CZKP4jg3XsVQh5QXNRrCmFXw1d9g+3uwbjoc/QFunwtVappdnYiIiEip8PHxwWq1cvz4cUJCQvDx8cFisZhdlgCGYZCWlsapU6ewWq34+Phc8bF0H53KaseHsPz/ICMZqtaGu+bBNdebXZWIiIhIqUhLS+PEiRMkJSWZXYrkwd/fn/Dw8DyDTmGzgYJOZRb3GyweAWf2gdULoqZB5wmgv2iIiIhIJWAYBhkZGTgcGspflthsNry8vPLtZVPQkcJJvQBfPAq/fOx63HQA3DoL/IJNLUtEREREJC+FzQaajKCys1eFO96Gm/8DVm/Ysxze7AEnfja7MhERERGRK6agI66hah3vhzErIegaOHcI3roRfpwHZb/DT0REREQkFwUdyVanHTywHq69CRypsHwifDoO0hLNrkxEREREpEgUdMSTf3UY+iFEPQMWG+z8COb2gVN7za5MRERERKTQFHQkN6sVbvg/GPkFBITCqd3wZi/YtdTsykRERERECkVBR/IX2RUe2ACR3SA9ET4eA8snQUaq2ZWJiIiIiBRIQUcKVjUURnwG3R53Pf7xbXinr2vCAhERERGRMkpBRy7PaoM+T8PdS8GvGhz/CeZ0h71fmV2ZiIiIiEieFHSk8Brf6BrKVqc9pMTDh0Nh1RRwZJhdmYiIiIiIBwUdKZrgCBj9FXQa73q86b/w7kBIOGFuXSIiIiIiOVxR0Jk1axaRkZH4+vrSqVMntmzZUmD7mTNn0qRJE/z8/IiIiOD//u//SElJuaKCpQzw8oF+L8Bd88GnKhz5DuZ0g4Prza5MRERERAS4gqCzaNEiJk2axNSpU9m+fTtt2rShb9++nDx5Ms/2Cxcu5Mknn2Tq1Kns3r2bt99+m0WLFvH3v//9qosXk7W4Dcaug9CWkHgK3h8E618Ep9PsykRERESkkrMYhmEUZYdOnTrRoUMHXnvtNQCcTicRERE8/PDDPPnkk7naP/TQQ+zevZvo6Gj3uscee4wffviBjRs3Fuo5ExISCAoKIj4+nsDAwKKUK6UhPRm+fBx+WuB63CgKbnsTqtQwty4RERERqXAKmw2K1KOTlpbGtm3biIqKyj6A1UpUVBSbN2/Oc58uXbqwbds29/C2gwcP8uWXX3LzzTfn+zypqakkJCR4LFKGefvBrbNci5cv7F/tGsp2dKvZlYmIiIhIJVWkoHP69GkcDgehoaEe60NDQ4mNjc1zn+HDh/Pss89yww034O3tTcOGDenZs2eBQ9emT59OUFCQe4mIiChKmWKW6/4C90VD9YaQcAzm3QTfvwFF6zQUEREREblqJT7r2rp163j++ed5/fXX2b59O5988gkrVqzgn//8Z777TJ48mfj4ePdy9OjRki5TiktYS9d1O80HgTMDvn4SFo9wTUctIiIiIlJKvIrSuGbNmthsNuLi4jzWx8XFERYWluc+Tz/9NPfccw/33XcfAK1atSIxMZGxY8fy1FNPYbXmzlp2ux273V6U0qQs8Q10zci25U1Y+RTs/hzifoHB70FYK7OrExEREZFKoEg9Oj4+PrRr185jYgGn00l0dDSdO3fOc5+kpKRcYcZmswFQxHkQpDyxWKDTA3Dv1xAUAWcPwltRsP09DWUTERERkRJX5KFrkyZNYu7cubz77rvs3r2b8ePHk5iYyOjRowEYMWIEkydPdrcfOHAgb7zxBh999BExMTGsWrWKp59+moEDB7oDj1RgddvDA99C4z9DRgp8/jAsexDSksyuTEREREQqsCINXQMYMmQIp06dYsqUKcTGxtK2bVu+/vpr9wQFR44c8ejB+cc//oHFYuEf//gHx44dIyQkhIEDB/Lcc88V36uQss2/OgxbBBtfhrXPwc8L4cTPMPhdqNnY7OpEREREpAIq8n10zKD76FQgMd/C0jGQeBJ8AuCW/0HL282uSkRERETKiRK5j47IVavfHcZtgHo3QNpFWDoavvwrZKSaXZmIiIiIVCAKOlL6qobBiM/ghkmux1vehHdugvNHzK1LRERERCoMBR0xh80Loqa6rt3xDYbj22F2N/h9pdmViYiIiEgFoKAj5mpyk2tWttp/gpTzsHAwrJ4GjgyzKxMRERGRckxBR8xXrZ7rfjsdx7oeb3wZ3rsVLsSaW5eIiIiIlFsKOlI2eNnh5hfhzndcs7Ed3ugayhazwezKRERERKQcUtCRsqXlHTB2HdRq7pqC+r1bYMNL4HSaXZmIiIiIlCMKOlL21GwM90VDm+FgOCH6WfhwCCSdNbsyERERESknFHSkbPLxh0Gvu24o6uUL+76BOd3hj21mVyYiIiIi5YCCjpRdFgv8aQTctxqqN4D4o/BOX/hhDhiG2dWJiIiISBmmoCNlX1gr13U7zW4BZzp89TdYMgpSEsyuTERERETKKAUdKR98g2Dwe9B3Oli94LdlMLcXxP1qdmUiIiIiUgYp6Ej5YbFA5wdh9FcQWAfO7Ie5feCnD8yuTERERETKGAUdKX8iOsIDG6BhH8hIhs8ehM8mQHqy2ZWJiIiISBmhoCPlU5UacPdS6PUPsFjhpwXwVhScOWB2ZSIiIiJSBijoSPlltUKPv8I9n0KVEIj7Beb0gF+XmV2ZiIiIiJhMQUfKvwY9XUPZrukCaRdgyUj46knISDO7MhERERExiYKOVAyB4TDyC+j6qOvxD2/A/Jvh/FFz6xIRERERUyjoSMVh84Ibn4WhH7qmo/5jK8zpBvtWmV2ZiIiIiJQyBR2peJreDA98C+FtIfkcfHAnrPkXOB1mVyYiIiIipURBRyqmapFw70poP8b1+NsX4f1BcPGkmVWJiIiISClR0JGKy9sXBrwMt78F3lUg5luY3Q0ObTK7MhEREREpYQo6UvG1vgvGroWQpnAxFt4dCBtfAafT7MpEREREpIQo6EjlENIE7l8DrYeC4YDVz8BHwyDprNmViYiIiEgJUNCRysOnCtw2Gwb+F2x2+P1r1w1Gj203uzIRERERKWYKOlK5WCzQbhSM+cY1YUH8EXinL2yZC4ZhdnUiIiIiUkwUdKRyqt0Wxq6HpgPAkQZfPg4fj4HUC2ZXJiIiIiLFQEFHKi+/YBiyAP78HFi94JeP4c1eEPeb2ZWJiIiIyFVS0JHKzWKBLg/BqBVQtTac2Qdze8PPH5ldmYiIiIhcBQUdEYBrrodxG6BBL8hIhk8fgM8fgfQUsysTERERkSugoCOSpUpN+MvH0HMyYIHt78LbUXD2oNmViYiIiEgRKeiI5GS1Qc8n4Z5PwL8GxO5yTUH92+dmVyYiIiIiRaCgI5KXhr1h3EaIuB5SE2DxPfD138GRbnZlIiIiIlIICjoi+QmsDaOWQ+eHXI+/nwXz+0P8MXPrEhEREZHLUtARKYjNG/o+B0M+AHsQHP0B5nSD/dFmVyYiIiIiBVDQESmMZgPggXUQ1hqSzsCCO2Dt8+B0mF2ZiIiIiORBQUeksKo3gDGroN1owID1M2DB7XDxlNmViYiIiMglFHREisLbFwbOhNveBG9/OLjONZTt8GazKxMRERGRHBR0RK5EmyFw/xqoeS1cOOGapGDTq2AYZlcmIiIiIijoiFy5Ws3g/rXQ6i4wHLDqafjobkg+b3ZlIiIiIpWego7I1bAHwO1zof/LYPOBvStgTnc4/pPZlYmIiIhUago6IlfLYoEOY2DMNxBcD84fhrf/DFvf1lA2EREREZMo6IgUl9rXwQProcnN4EiDFZPgk7GQetHsykREREQqHQUdkeLkVw2GLoQbnwWLDXYthrm94eQesysTERERqVQUdESKm8UCXR+FUcshIAxO74W5vWDnYrMrExEREak0rijozJo1i8jISHx9fenUqRNbtmwpsP358+eZMGEC4eHh2O12rr32Wr788ssrKlik3KjXBcZthPo9ID0JPrkfvpgI6SlmVyYiIiJS4RU56CxatIhJkyYxdepUtm/fTps2bejbty8nT57Ms31aWho33ngjhw4dYunSpezdu5e5c+dSp06dqy5epMwLCIF7PoXufwMssG0evPNnOBtjdmUiIiIiFZrFMIo2LVSnTp3o0KEDr732GgBOp5OIiAgefvhhnnzyyVztZ8+ezYsvvsiePXvw9va+oiITEhIICgoiPj6ewMDAKzqGiOn2rXb16iSfBXsQ3PYGNO1vdlUiIiIi5Uphs0GRenTS0tLYtm0bUVFR2QewWomKimLz5s157vP555/TuXNnJkyYQGhoKC1btuT555/H4XDk+zypqakkJCR4LCLlXuMoGLcB6naE1Hj4aDh88w9wpJtdmYiIiEiFU6Sgc/r0aRwOB6GhoR7rQ0NDiY2NzXOfgwcPsnTpUhwOB19++SVPP/00L730Ev/617/yfZ7p06cTFBTkXiIiIopSpkjZFVQXRq2A6ye4Hn/3P3h3ICQcN7cuERERkQqmxGddczqd1KpVizfffJN27doxZMgQnnrqKWbPnp3vPpMnTyY+Pt69HD16tKTLFCk9Xj5w0/Mw+D2wB8KRzTC7GxxYa3ZlIiIiIhVGkYJOzZo1sdlsxMXFeayPi4sjLCwsz33Cw8O59tprsdls7nXNmjUjNjaWtLS0PPex2+0EBgZ6LCIVTvNbYew6CG0FSafh/dtg3QxwOs2uTERERKTcK1LQ8fHxoV27dkRHR7vXOZ1OoqOj6dy5c577dO3alf379+PM8eHt999/Jzw8HB8fnyssW6SCqNEQ7lsFfxoBGLDuefjgDkg8bXZlIiIiIuVakYeuTZo0iblz5/Luu++ye/duxo8fT2JiIqNHjwZgxIgRTJ482d1+/PjxnD17lkcffZTff/+dFStW8PzzzzNhwoTiexUi5Zm3H9zyPxj0Bnj5wYE1rqFsR34wuzIRERGRcsurqDsMGTKEU6dOMWXKFGJjY2nbti1ff/21e4KCI0eOYLVm56eIiAhWrlzJ//3f/9G6dWvq1KnDo48+yhNPPFF8r0KkImg7HMLbwOKRcGYfzL8ZoqZB5wlgsZhdnYiIiEi5UuT76JhB99GRSiX1Anz+CPz6ietx0wEw6HXwDTK3LhEREZEyoETuoyMipcBeFe58B27+D1i9Yc9ymNMDTuw0uzIRERGRckNBR6Qsslig4/1w70oIugbOxcBbUbBtPpT9TlgRERER0ynoiJRlddvBA+uhcV9wpMIXj8Kn4yAt0ezKRERERMo0BR2Rss6/Ogz7CPpMBYsVdn4Ec/vAqd/NrkxERESkzFLQESkPrFboNglGfgEBoXBqN7zZE3YtNbsyERERkTJJQUekPIm8AR7YAJHdID0RPh4DKx6DjFSzKxMREREpUxR0RMqbqqFwzzLo9pjr8da34J2+cO6wqWWJiIiIlCUKOiLlkc0L+kyB4UvArxoc/wnmdIO9X5ldmYiIiEiZoKAjUp5d+2fXULY67SElHj4cCqumgiPD7MpERERETKWgI1LeBUfA6K+g0zjX400z4b1b4EKsqWWJiIiImElBR6Qi8PKBfjPgrvngUxUOb4LZN8DB9WZXJiIiImIKBR2RiqTFbTB2HdRqAYmn4P1BsP5FcDrNrkxERESkVCnoiFQ0NRvBfauh7V/AcMLaf8HCuyDxjNmViYiIiJQaBR2RisjHHwbNgltngZcv7F8Nc7rD0a1mVyYiIiJSKhR0RCqy6/7i6t2p3hAS/oB5/eD7N8AwzK5MREREpEQp6IhUdGGtXNftNL8VnOnw9ZOwZCSkJJhdmYiIiEiJUdARqQx8A+Gud+GmGWD1ht8+gzd7QOwusysTERERKREKOiKVhcUC14+De7+GwLpw9iC8FQXb3ze7MhEREZFip6AjUtnUbQ/jNkCjGyEjBT5/CJY9CGlJZlcmIiIiUmwUdEQqI//qMHwx9P4HWKyw4wNX787p/WZXJiIiIlIsFHREKiurFbr/FUZ8BlVqwclfXdft/PKJ2ZWJiIiIXDUFHZHKrn5311C2el0h7SIsHQ1f/g0y0syuTEREROSKKeiICFQNgxGfww3/53q8ZQ7MuwnOHzG3LhEREZErpKAjIi42L4h6BoYtAt9gOLYN5nSH378xuzIRERGRIlPQERFPTW6CB76F2tdB8jlYeBesngaODLMrExERESk0BR0Rya1aPbh3JXS43/V448vw/iC4EGdqWSIiIiKFpaAjInnzskP//8Cd74BPABzaAHO6waGNZlcmIiIiclkKOiJSsJZ3wNh1ENIMLsbBuwNhw0vgdJpdmYiIiEi+FHRE5PJqNob7o6HNMDCcEP0sfDgUks6aXZmIiIhInhR0RKRwfKrAoDdg4Ktgs8O+lTCnB/yxzezKRERERHJR0BGRwrNYoN1IuG81VKsP8Ufgnb7ww5tgGGZXJyIiIuKmoCMiRRfeGh5YD80GgjMdvvorLB0NqRfMrkxEREQEUNARkSvlGwSD34e+z4PVC379FN7sCXG/ml2ZiIiIiIKOiFwFiwU6T4BRX0JgHTizH+b2gR0Lza5MREREKjkFHRG5etd0ggc2QMM+kJEMy8bDZw9BerLZlYmIiEglpaAjIsWjSg24eyn0egqwwE/vw1s3wpkDZlcmIiIilZCCjogUH6sVevwN7vkU/GtC3C7XFNS/fWZ2ZSIiIlLJKOiISPFr2AvGbYBrOkPaBVg8Ar6eDBlpZlcmIiIilYSCjoiUjMDaMPIL6PKI6/H3r8P8myH+D3PrEhERkUpBQUdESo7NG/78Txj6oWs66j+2wuxusG+12ZWJiIhIBaegIyIlr+nN8MC3EN4Gks/CB3fCmn+B02F2ZSIiIlJBKeiISOmoFgn3fgPt7wUM+PZFeP82uHjS7MpERESkAlLQEZHS4+0LA16B298C7yoQs941lO3wd2ZXJiIiIhWMgo6IlL7Wd8HYtRDSFC7GwvwBsHEmOJ1mVyYiIiIVhIKOiJgjpAncvwZaDwHDAaunwkfDIfmc2ZWJiIhIBXBFQWfWrFlERkbi6+tLp06d2LJlS6H2++ijj7BYLAwaNOhKnlZEKhqfKnDbHBgwE2w+8PtXMKc7HP/J7MpERESknCty0Fm0aBGTJk1i6tSpbN++nTZt2tC3b19Oniz4guJDhw7x+OOP061btysuVkQqIIsF2o+GMasguB6cPwJv/xm2vgWGYXZ1IiIiUk4VOei8/PLL3H///YwePZrmzZsze/Zs/P39eeedd/Ldx+FwcPfddzNt2jQaNGhwVQWLSAVVu61rCuqmA8CRBiseg4/vg9SLZlcmIiIi5VCRgk5aWhrbtm0jKioq+wBWK1FRUWzevDnf/Z599llq1arFmDFjCvU8qampJCQkeCwiUgn4BcOQBfDnf4HFBr8shbm94ORusysTERGRcqZIQef06dM4HA5CQ0M91oeGhhIbG5vnPhs3buTtt99m7ty5hX6e6dOnExQU5F4iIiKKUqaIlGcWC3R5GEatgKrhcPp3mNsbfl5kdmUiIiJSjpTorGsXLlzgnnvuYe7cudSsWbPQ+02ePJn4+Hj3cvTo0RKsUkTKpHqd4YEN0KAnpCfBp2Phi0chPcXsykRERKQc8CpK45o1a2Kz2YiLi/NYHxcXR1hYWK72Bw4c4NChQwwcONC9zpl5nwwvLy/27t1Lw4YNc+1nt9ux2+1FKU1EKqKAEPjLJ7D+37B+BmybD8e2w+B3obqu9xMREZH8FalHx8fHh3bt2hEdHe1e53Q6iY6OpnPnzrnaN23alF27drFjxw73csstt9CrVy927NihIWkicnlWG/SaDH/5GPxrQOxOmNMTdn9hdmUiIiJShhWpRwdg0qRJjBw5kvbt29OxY0dmzpxJYmIio0ePBmDEiBHUqVOH6dOn4+vrS8uWLT32Dw4OBsi1XkSkQI36uIayLR0NR3+ARX+Bzg9B1DNg8za7OhERESljihx0hgwZwqlTp5gyZQqxsbG0bduWr7/+2j1BwZEjR7BaS/TSHxGprILquCYpWP0MbH7NtfzxI9z5jmubiIiISCaLYZT9O/IlJCQQFBREfHw8gYGBZpcjImXB7i9g2YOQmuAa0nbHW9Cwt9lViYiISAkrbDZQ14uIlE/NBsID6yGsNSSdgfdvh7XTwekwuzIREREpAxR0RKT8qt4AxqyCdqMAA9a/AAtuh8TTZlcmIiIiJlPQEZHyzdsXBv4XbpsD3v5wcB3M7gZHvje7MhERETGRgo6IVAxthsL9a6DmtXDhOMy7GVY8DvtXQ3qy2dWJiIhIKdNkBCJSsaRehC8ehV+WZq/z8oPIG6BRlGup0RAsFvNqFBERkStW2GygoCMiFY9hwL5vYM9y2B8NCcc8twfXyw499buDPcCcOkVERKTIFHRERMAVek7tcQ1h27cKjmwGR1r2dqs31OucHXxqNVdvj4iISBmmoCMikpfUi3Booyv47F8F5w55bq9aGxr1cYWeBj3BL9iEIkVERCQ/CjoiIoVx5kBm6FkNMRsgI8fEBRYb1O3gCj2NoyCsDVg1h4uIiIiZFHRERIoqPQWOfAf7MoPP6b2e2/1rZvf2NOwNVWqaU6eIiEglpqAjInK1zh9xTWawfzUcXA9pF3JstEDt67Kv7anTDmxeppUqIiJSWSjoiIgUp4w0+GNL9jC32F2e232DoWGvzN6ePhAYbkqZIiIiFZ2CjohISboQCwfWuGZyO7AGUs57bg9tmT3MLeJ68PIxpUwREZGKRkFHRKS0OB1wbHv2TG7HtgM5/mv1CYD6PbKDT7V6ppUqIiJS3inoiIiYJfEMHFybPcwt8ZTn9hqNofGNruBTryt4+5lTp4iISDmkoCMiUhY4nRC3yzXEbX80HP0BDEf2di9fiLwhc1KDG6FGQ92wVEREpAAKOiIiZVFKvGsGt6zenoRjntuD62XP5Fa/O9gDzKlTRESkjFLQEREp6wwDTu3JDj2HvwNHWvZ2qzfU65wdfGo1V2+PiIhUego6IiLlTepFOLQxe1KDc4c8t1cNz5zQ4EZo0BP8gk0oUkRExFwKOiIi5d2ZA5k3LF0FMRsgIzl7m8UGdTtk9vb0gfC2YLWaVqqIiEhpUdAREalI0lPgyHeZwWe1a8hbTv41s6evbtgbqtQ0p04REZESpqAjIlKRnT+SHXoOroe0Czk2WqB22+yZ3Oq0A5uXWZWKiIgUKwUdEZHKwpEOR7e4hrjtXw2xuzy3+wZBg17ZkxoEhptTp4iISDFQ0BERqawuxMKBNa7Qc2ANJJ/z3B7aMnuYW8T14OVjTp0iIiJXQEFHRETA6YBj27OnsD62Dcjx375PgOt+PVm9PdXqmVaqiIhIYSjoiIhIbklns3t79q+GxFOe22s0zg49kV3B28+cOkVERPKhoCMiIgVzOiFulyvw7FsNR38Aw5G93csXIm/IDj41GumGpSIiYjoFHRERKZqUeNcMblm9PQnHPLcH18sOPfW7gb2qOXWKiEilpqAjIiJXzjDg1N7smdwOfweOtOztVm+45npX6Gl8I9Rqrt4eEREpFQo6IiJSfNIS4dDGzGFuq+BcjOf2quHZM7k16AV+waaUKSIiFZ+CjoiIlJwzB7JvWBrzLWQkZ2+z2KBuh8xhbn0gvC1YraaVKiIiFYuCjoiIlI70FDjyXXbwObXHc7t/TWjY2zXErWFvqFLTnDpFRKRCUNARERFznD8KB6JdQ9wOroe0Czk2WqB22+xJDeq0B5uXWZWKiEg5pKAjIiLmc6TD0S3ZM7nF7vTc7hvkuqYna5hbYG1z6hQRkXJDQUdERMqeC3Gu3p79q103Lk0+57m9VgtonNnbE3E9ePmYU6eIiJRZCjoiIlK2OR1w/CfXELf9q+HYNiDHrySfAKjfPXs2t2qRZlUqIiJliIKOiIiUL0lnXb08WZMaJJ703F6jcfa1PZFdwdvPnDpFRMRUCjoiIlJ+OZ0Qtyvz2p5oOPI9GI7s7V6+UK+raya3RlFQo5FuWCoiUkko6IiISMWREu+6X8++Va7gk/CH5/bgazJ7e26E+t3AXtWcOkVEpMQp6IiISMVkGHBqb2Zvzyo4/B040rK3W73hmuuzh7mFtlBvj4hIBaKgIyIilUNaIhzamD2F9dmDnturhmdPaNCgJ/hVM6VMEREpHgo6IiJSOZ05kD2hwaENkJ6Uvc1ihbodsnt7wtuC1WpaqSIiUnQKOiIiIukpcGRzdm/PqT2e2/1rQMM+2TcsrVLTnDpFRKTQFHREREQudf5o9g1LD66H1IQcGy1Qu212b0+d9mDzMqtSERHJh4KOiIhIQRzpcHRLdm9P7E7P7fYgaNjTNZNboz4QWNuUMkVExFNhs8EVDUyeNWsWkZGR+Pr60qlTJ7Zs2ZJv27lz59KtWzeqVatGtWrViIqKKrC9iIhIqbB5u248GjUVxm2Ax36HQbOh5Z2uCQtS4+G3z+Dzh+DlZvB6F/jmaVdPUEba5Y8vIiKmKnKPzqJFixgxYgSzZ8+mU6dOzJw5kyVLlrB3715q1aqVq/3dd99N165d6dKlC76+vsyYMYNPP/2UX3/9lTp16hTqOdWjIyIipcrpgOM/Zff2/PEjkOPXpXcVaNAjeza3apFmVSoiUumU2NC1Tp060aFDB1577TUAnE4nERERPPzwwzz55JOX3d/hcFCtWjVee+01RowYUajnVNARERFTJZ2FA2uyZ3NLPOm5vUajzCFuUa5eIm8/c+oUEakECpsNinSVZVpaGtu2bWPy5MnudVarlaioKDZv3lyoYyQlJZGenk716tXzbZOamkpqaqr7cUJCQr5tRURESpx/dWh1p2txOiHuF9fNSvdHw9Ef4Mx+1/LDG+DlC/W6ukJP4xtdIUg3LBURKXVFCjqnT5/G4XAQGhrqsT40NJQ9e/bks5enJ554gtq1axMVFZVvm+nTpzNt2rSilCYiIlI6rFYIb+1auj0GKfEQ862rp2ffakj4wzWz24FoWDkZgq/JnsmtfnewVzX7FYiIVAqlOm/mCy+8wEcffcS6devw9fXNt93kyZOZNGmS+3FCQgIRERGlUaKIiEjR+AZBs4GuxTDg1N7sa3sOb4LzR+DHd1yL1RuuuT47+IS2UG+PiEgJKVLQqVmzJjabjbi4OI/1cXFxhIWFFbjvf/7zH1544QVWr15N69atC2xrt9ux2+1FKU1ERMR8FgvUaupaujwEaYlwaGN28Dl7EA5tcC2rp0LVcNcNSxtHQYOertneRESkWBQp6Pj4+NCuXTuio6MZNGgQ4JqMIDo6moceeijf/f7973/z3HPPsXLlStq3b39VBYuIiJQbPlXg2r6uBeDMgcxJDVa7hrtdOAE7FrgWixXqdsjs7ekD4de5hsmJiMgVuaLppUeOHMmcOXPo2LEjM2fOZPHixezZs4fQ0FBGjBhBnTp1mD59OgAzZsxgypQpLFy4kK5du7qPExAQQEBAQKGeU7OuiYhIhZOeAkc2Z/b2RMOp3Z7b/Wu4ensaRUHD3hAQYk6dIiJlTIlNLw3w2muv8eKLLxIbG0vbtm159dVX6dSpEwA9e/YkMjKS+fPnAxAZGcnhw4dzHWPq1Kk888wzxfpiREREyq34P7KHuB1cD6mXzDga3tY1i1ujKKjTHmylepmtiEiZUaJBp7Qp6IiISKXiSIc/tsK+Va7gE7vTc7s9CBr2zOzt6QNBhbsBt4hIRaCgIyIiUlFciMu+tudANCSf89xeq4Xrup5GUa5Z3bw0oY+IVFwKOiIiIhWR0wHHf8oe5vbHj0COX+XeVaBBj+zgUy3SrEpFREqEgo6IiEhlkHQWDq513ax0/2pIPOm5vUajzJncboTIruDtZ06dIiLFREFHRESksnE6Ie6X7Jncjn4Pzozs7V6+UK9r9g1LazbWDUtFpNxR0BEREansUuJd9+vZv9rV45Pwh+f24GuyQ0/97mCvak6dIiJFoKAjIiIi2QwDTv+eGXpWweFN4EjL3m71dk1kkBV8Qluot0dEyiQFHREREclfWiIc2pQ5zG0VnD3ouT0gLDP09IGGvcCvmjl1iohcQkFHRERECu/MgewprGO+hfSk7G0WK9TtkB18wq8Dq9W8WkXEFIZhYCkDPb0KOiIiInJlMlLh8HfZkxqc2u253b8GNOztmsmtYW8ICDGnTpES4nQaOAwDh9O1ZDizv3Y9duJ0QobT6VpnGGQ4DPfXDqfrsdPI2teJwwkOp9PjWBlOA2eex8/c15G5b9ZxHJ515XcM93Mamc/pUUsh9720FqfBA90b8mS/pmafnkJnA69SrElERETKAy+7a7haw17Q9zmI/8MVePavhoPrIOkM7FriWgDC22Zf21O3A9j08aK8KFMf6A3DvV++bUvpA33Z7wYwh7OcfWPUoyMiIiKF50iHP7Zm37D0xM+e2+1B0LCnK/Q07ANBdUwp81Il+YHekXO/PP8i73kMMz7Qu2pxetaiD/RXxMtqwZZj8XhssWCzWfCyWrFawMtqzdXWesk+nsezutpYMtfbMo+Zxz7WXMezYrOAzeY6hns/W47jeRzDitXqWWNer81qcR3DZrHgb/ciwG7+HzI0dE1ERERK3oU4nPujcexbje3gWqwpZz02Xwy6lthaN/BH9a4crdqaixk2ktMySE53kJTmIDndkR0GLvOBPr9QoA/0JaOgD/S5PiRf9kP1pftaivyB3ma1YrOSHQZyts36MJ5nLUX8QJ91PKs1V11Wq/nXp4iCjoiIiGRyOI3MYJFBSpqTpPQMktIcpKS5wkZSetbXGTm+dnh8nZzuIDnNQVJ6BslpWV+7/k3NcAJgxUkry0F6WHfS07aDNpYD2CzZHzMSDTvfOVuw3tmGdc42/GHUMutb4sGrgA/Vl/5F/ko/0Hv+tT/3B/yCPsgX5QN9/rXk7knIqjGrLn2gl/JCQUdERKSccDgNkjJ7OZIvDRburzM81ucMGtlf5+gpyWyblOYgLTOIlDSLBfy8ba7Fx0aoVyLXs4sOGdtpm7aNYIdnb89Zv2s4Vr0LCVUbkuobQqpfLVJ9Q0j3q4nNyyefv7ZbCwwbhflAnzNo6AO9SPmjyQhERESKSYbDmSt4eIaJDFIy1yWlOfL4OsNjfc4Qkpxe+kHE38eGb+a/fj5e+Hlb8ffxws/H5t7u+bVX9teZIcbza9d2X29rHlPPDnT9YxgQuyt7Jrej31M9+QjVjx3Jq1KoUtN1L5+qYVA1NMfXYVA1HAJCXYuXT0l/20SknFKPjoiIlHsZDqdn74ZHj0gRe0ryaJvmKL0g4p8ZHvx8bPh7e+HrY8M/K5zk+toLPx8rfj5el+yXI4zkCCl2r7yCiElSElz364n51jWr28VYuBALF+PAmVH44/jXyAxBodkBqGq4ZzgKCAVv35J7LSJSqjR0TUREyoz0S3tE0hwkZ14nkpxPD0lympPkzOtB8g4o2SEm3VE6v8qyg4hXvr0befaC5NlT4pUrmJSpIGIWp9M1ffXFWLgQBxdOZIegrCCU9bUzvfDH9Q3ODkDuQJQVhHL0Gvn4l9hLE5HioaFrIiJSaOkOZx5BI68AkvdQrMsN1yqtIGK1kGsIVtYQLc+vvfJZn/W1V54hRkGkFFitrhuQBoRAWKv82xkGJJ9zBSF3CMoZjnIEIkcqpJx3LZfe/PRS9qB8hsuFeT72qVKcr1pESoCCjohIOZCW4fQIIR7Xe1xmKFbuHhBnrovWM5ylH0Ty7hHJewhWYXtKfGwKIpWGxQL+1V1LaIv82xmGK+B49AidcAWinD1FF2IhIxlS413L6b0FP79P1QJ6h3IMmbNXddUqIqVOQUdE5CoZhkGaw+metjfvoVauoVhF7SnJ2l5aQcRmteQTLmz4eXvlCh8eX+dsc0mIcV1rYlUQkdJnsYBfNddSq1n+7QwDUhPy7hG6NBClJ0LaBThzAc7sL/j5vasUPFwua71vkAKRSDFT0BGRCi8riOQXMHJNyZvntL15X1OSkrndUUpBxMtqyWdY1uWHYnnOtOUZYvwzA4q3zaIgIpWTxeIKG75BEHJtwW1TL+TRI3RpOIpzBaf0RDh70LUUxMsv93C5nBMrZAUiv2oKRCKFpKAjIuVWcpqD4/HJHD+fzInzKRw7n8yJ+GSOn0/heHwyF1Kye09KO4h4DMXKd0re7DaFumjd24aPl7VUXoeIFMBe1bXUbFRwu7TEgofLZa1PiXcNmzt3yLUUxGbPEYjym2ku3DWkT4FIKjkFHREpkzIcTk5eSOX4+WSOx6dkhplkjp3P/Do+mXNJRZhxKZO3zZLr3h95BZAr7SnxtimIiEgmnypQo6FrKUh6cu4Z5fIKRMnnXBMrnD/iWgpi9c7RK5TPcLmq4a7pua36f0sqJgUdESl1hmFwPik9swcmJTPMZPbEZAaauAupheqFqeJjo3awX+biS+0g19fhwb4E+/nkuJeIK8goiIhImePtB9Xru5aCZKR6hqH8ZppLOu2aejv+qGspiNULqtQqeLhc1TCoEgJWW/G9ZpFSoKAjIsUua0jZiczgcumQsuPnk0lJv/wNGL2sFsKCfKkd7EedYD/CM7+uHez6NzzIj0BfL11TIiKVg5cdgq9xLQXJSIPEkwX3Dl2Ig8RTrpuzXjjuWgpisWYGogJuzFo1zNXGpo+XUjbonSgiRVLQkLITmSGmsEPKagb4uIJLkKsHpk5meMkKMjUD7NisCjEiIkXi5QNBdV1LQRwZeQSinL1DWYHoJBhO1/aLsXDi5wIOanH1/uQ701xWOAoFm3exvmyRSynoiIhbQUPKTpx3hZirGVIWnuPrsCBffL01DEJExDQ2Lwis7VoK4nS4en8u7RHKa6Y5w+EKRoknIXZXwcf1r1nwcLmsbV724nvNUqko6IhUIpcOKXPPWBafOWPZ+RSS0x2XPU7OIWW1M/8ND/ajTrBvZo+MhpSJiFQYVlt28CiI0wlJZ3L3COUaOhfruoYo6bRriful4OP6VSt4uFxWMPL2K77XLBWCgo5IBZE1pOxEfOYwssweGPfwsvgUziamFepYWUPK3NfEBPl5XBujIWUiIpKL1QoBIa6lIE6nawa5Cyfynkwh53VFjjRX2+RzcPK3go/rG3T5G7NWDXPNhieVgoKOSDmQNaTMPYwsPtndA3P8CoeUZfXAaEiZiIiUKqsVqtRwLbTMv51huAJOfr1DOa8rykh23Y8oJR5O7Sn4+e2BuYfHVQ3PHY7sVYv1ZUvpU9ARKQNS0h2ZgSV7VrKrGlKW44J+DSkTEZFyyWJx3fjUvzrUapZ/O8NwBZzL3Zj1QhykJ0Jqgms5s6/g5/cJuMxwucz19kDdnLWMUtARKWEOp8HJC1nTLF/9kLKcs5LlvGdMHQ0pExGRyshiAb9g1xLSpOC2qRfy6BHKIxClXYC0i3D2Ipw9UPAxvf0vf2PWqqHgG6xAVMoUdESuwuWGlJ2ITyE2IeWKhpRl9cBkXeyvIWUiIiJXyV7VtdRsXHC71IvZ1wzlN1zuQiykxkN6EpyLcS0F8fLNDkT5zjQX7pp8QYGoWCjoiBQgryFlJzyGl13ZkLLwYM8QUzvIj0A/DSkTEREpE+wBrqVGw4LbpSXlHYguHUaXfA4yUuD8YddSEJtPdq9QrqFzOcKRX3XX9U6SLwUdqbRyDik7nqMHxnUPGde6KxlSFh7k57rxZY7hZSFVNaRMRESkwvHxh+r1XUtB0lMuCUT5zDSXdMY101z8EddSEKtXjh6ifIbLBYRBlZquKcIrIQUdqZAMwyA+OT17GFl8dpgp6pAy/5w3vsy6Z0yQ65oYDSkTERGRy/L2hWr1XEtBMtJc4cfdI3TpTVozA1HiaXBmQMIx11IQiw0Cal3mxqxhUCXEdRPZCqRivRqpNLKGlGX1wFzNkLLQQF/PHhgNKRMREREzePlAcIRrKYgjHS6eLHi43IVYSDwFhiNz/Qk4sSP/Y1qsrrCT70xzYa6gVqVmsb7kkqSgI2VOcQ4pq1HFx32jSw0pExERkQrB5g1BdVxLQRwZkHQ6R49QPjdmvXjSFYiyepNid+Z9vI5j4eYXi//1lBAFHSlVWUPK3MPI4nOGGdfXxTGkLDzzaw0pExERkUrL5pXdI1MQp8N1fdDlbswaVLd06i4mCjpSrC4dUnbCI9C41ielXX5Imc1qIezSIWXuMOPqmdGQMhEREZFiYM26jqcWhJtdTPFR0JFCyx5S5tkDkxVkTpxP4UwRh5SFZ10LE+wZYjSkTERERESuhoKOALmHlLlufJmSGWZcgSYuIYWMIgwpcw8jC8oOMrU1pExERERESoGCTiWRku7gRHz29MrHr3JImWcPjIaUiYiIiEjZoqBTATicBqcupOaYlezqhpSFB/tSO8hPQ8pEREREpNxS0CnjDMMgITkj+14xVzmkLOuaGA0pExEREZGK7IqCzqxZs3jxxReJjY2lTZs2/O9//6Njx475tl+yZAlPP/00hw4donHjxsyYMYObb775iouuSPIaUuYKM8nu9VczpCw8R89MkJ+3hpSJiIiISKVQ5KCzaNEiJk2axOzZs+nUqRMzZ86kb9++7N27l1q1auVq/9133zFs2DCmT5/OgAEDWLhwIYMGDWL79u20bNmyWF5EWVXQkLKsEHM1Q8pyhphaVX01pExEREREJJPFMIzLj3nKoVOnTnTo0IHXXnsNAKfTSUREBA8//DBPPvlkrvZDhgwhMTGR5cuXu9ddf/31tG3bltmzZxfqORMSEggKCiI+Pp7AwMCilFticg4pc4eYS3pmCjukzM/blj2ELDO8hAdnzVjmWq8hZSIiIiIihc8GRerRSUtLY9u2bUyePNm9zmq1EhUVxebNm/PcZ/PmzUyaNMljXd++fVm2bFm+z5Oamkpqaqr7cUJCQlHKLDHHzyfzxMc7r3hIWc4emNpB2WFGQ8pERERERIpXkYLO6dOncTgchIaGeqwPDQ1lz549ee4TGxubZ/vY2Nh8n2f69OlMmzatKKWVCj9vGxv2nfZYV72KjzvE1PG4AaaGlImIiIiImKVMzro2efJkj16ghIQEIiIiTKzIJdjfmxfvbO2eoSw8yA8/Hw0pExEREREpa4oUdGrWrInNZiMuLs5jfVxcHGFhYXnuExYWVqT2AHa7HbvdXpTSSoXFYuGu9uYHLhERERERKZi1KI19fHxo164d0dHR7nVOp5Po6Gg6d+6c5z6dO3f2aA+watWqfNuLiIiIiIhcrSIPXZs0aRIjR46kffv2dOzYkZkzZ5KYmMjo0aMBGDFiBHXq1GH69OkAPProo/To0YOXXnqJ/v3789FHH/Hjjz/y5ptvFu8rERERERERyVTkoDNkyBBOnTrFlClTiI2NpW3btnz99dfuCQeOHDmC1ZrdUdSlSxcWLlzIP/7xD/7+97/TuHFjli1bVuHvoSMiIiIiIuYp8n10zFAW76MjIiIiIiKlr7DZoEjX6IiIiIiIiJQHCjoiIiIiIlLhKOiIiIiIiEiFo6AjIiIiIiIVjoKOiIiIiIhUOAo6IiIiIiJS4RT5PjpmyJoBOyEhweRKRERERETETFmZ4HJ3ySkXQefChQsAREREmFyJiIiIiIiUBRcuXCAoKCjf7eXihqFOp5Pjx49TtWpVLBaLqbUkJCQQERHB0aNHdfPSCkLntGLSea14dE4rJp3XikfntOIpa+fUMAwuXLhA7dq1sVrzvxKnXPToWK1W6tata3YZHgIDA8vEiZbio3NaMem8Vjw6pxWTzmvFo3Na8ZSlc1pQT04WTUYgIiIiIiIVjoKOiIiIiIhUOAo6RWS325k6dSp2u93sUqSY6JxWTDqvFY/OacWk81rx6JxWPOX1nJaLyQhERERERESKQj06IiIiIiJS4SjoiIiIiIhIhaOgIyIiIiIiFY6CjoiIiIiIVDgKOnmYNWsWkZGR+Pr60qlTJ7Zs2VJg+yVLltC0aVN8fX1p1aoVX375ZSlVKoVVlHM6f/58LBaLx+Lr61uK1crlfPvttwwcOJDatWtjsVhYtmzZZfdZt24df/rTn7Db7TRq1Ij58+eXeJ1SNEU9r+vWrcv1s2qxWIiNjS2dguWypk+fTocOHahatSq1atVi0KBB7N2797L76fdq2XUl51S/V8u+N954g9atW7tvCNq5c2e++uqrAvcpDz+nCjqXWLRoEZMmTWLq1Kls376dNm3a0LdvX06ePJln+++++45hw4YxZswYfvrpJwYNGsSgQYP45ZdfSrlyyU9Rzym47vx74sQJ93L48OFSrFguJzExkTZt2jBr1qxCtY+JiaF///706tWLHTt2MHHiRO677z5WrlxZwpVKURT1vGbZu3evx89rrVq1SqhCKar169czYcIEvv/+e1atWkV6ejp//vOfSUxMzHcf/V4t267knIJ+r5Z1devW5YUXXmDbtm38+OOP9O7dm1tvvZVff/01z/bl5ufUEA8dO3Y0JkyY4H7scDiM2rVrG9OnT8+z/eDBg43+/ft7rOvUqZPxwAMPlGidUnhFPafz5s0zgoKCSqk6uVqA8emnnxbY5m9/+5vRokULj3VDhgwx+vbtW4KVydUozHldu3atARjnzp0rlZrk6p08edIAjPXr1+fbRr9Xy5fCnFP9Xi2fqlWrZrz11lt5bisvP6fq0ckhLS2Nbdu2ERUV5V5ntVqJiopi8+bNee6zefNmj/YAffv2zbe9lK4rOacAFy9epF69ekRERBT4Fw0pH/RzWrG1bduW8PBwbrzxRjZt2mR2OVKA+Ph4AKpXr55vG/28li+FOaeg36vlicPh4KOPPiIxMZHOnTvn2aa8/Jwq6ORw+vRpHA4HoaGhHutDQ0PzHfMdGxtbpPZSuq7knDZp0oR33nmHzz77jAULFuB0OunSpQt//PFHaZQsJSC/n9OEhASSk5NNqkquVnh4OLNnz+bjjz/m448/JiIigp49e7J9+3azS5M8OJ1OJk6cSNeuXWnZsmW+7fR7tfwo7DnV79XyYdeuXQQEBGC32xk3bhyffvopzZs3z7Ntefk59TK7AJGypnPnzh5/wejSpQvNmjVjzpw5/POf/zSxMhHJqUmTJjRp0sT9uEuXLhw4cIBXXnmF999/38TKJC8TJkzgl19+YePGjWaXIsWksOdUv1fLhyZNmrBjxw7i4+NZunQpI0eOZP369fmGnfJAPTo51KxZE5vNRlxcnMf6uLg4wsLC8twnLCysSO2ldF3JOb2Ut7c31113Hfv37y+JEqUU5PdzGhgYiJ+fn0lVSUno2LGjflbLoIceeojly5ezdu1a6tatW2Bb/V4tH4pyTi+l36tlk4+PD40aNaJdu3ZMnz6dNm3a8N///jfPtuXl51RBJwcfHx/atWtHdHS0e53T6SQ6OjrfMYqdO3f2aA+watWqfNtL6bqSc3oph8PBrl27CA8PL6kypYTp57Ty2LFjh35WyxDDMHjooYf49NNPWbNmDfXr17/sPvp5Lduu5JxeSr9Xywen00lqamqe28rNz6nZsyGUNR999JFht9uN+fPnG7/99psxduxYIzg42IiNjTUMwzDuuece48knn3S337Rpk+Hl5WX85z//MXbv3m1MnTrV8Pb2Nnbt2mXWS5BLFPWcTps2zVi5cqVx4MABY9u2bcbQoUMNX19f49dffzXrJcglLly4YPz000/GTz/9ZADGyy+/bPz000/G4cOHDcMwjCeffNK455573O0PHjxo+Pv7G3/961+N3bt3G7NmzTJsNpvx9ddfm/USJA9FPa+vvPKKsWzZMmPfvn3Grl27jEcffdSwWq3G6tWrzXoJconx48cbQUFBxrp164wTJ064l6SkJHcb/V4tX67knOr3atn35JNPGuvXrzdiYmKMnTt3Gk8++aRhsViMb775xjCM8vtzqqCTh//973/GNddcY/j4+BgdO3Y0vv/+e/e2Hj16GCNHjvRov3jxYuPaa681fHx8jBYtWhgrVqwo5YrlcopyTidOnOhuGxoaatx8883G9u3bTaha8pM1rfClS9Z5HDlypNGjR49c+7Rt29bw8fExGjRoYMybN6/U65aCFfW8zpgxw2jYsKHh6+trVK9e3ejZs6exZs0ac4qXPOV1PgGPnz/9Xi1fruSc6vdq2Xfvvfca9erVM3x8fIyQkBCjT58+7pBjGOX359RiGIZRev1HIiIiIiIiJU/X6IiIiIiISIWjoCMiIiIiIhWOgo6IiIiIiFQ4CjoiIiIiIlLhKOiIiIiIiEiFo6AjIiIiIiIVjoKOiIiIiIhUOAo6IiIiIiJS4SjoiIiIiIhIhaOgIyIiIiIiFY6CjoiIiIiIVDgKOiIiIiIiUuH8P91INOocEiE5AAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 1000x400 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Solve the factor graph and retrieve the results\n",
        "result = graph.optimize()\n",
        "\n",
        "states = []\n",
        "controls = []\n",
        "for k in range(N):\n",
        "    xk = X(k)\n",
        "    uk = U(k)\n",
        "    states.append(result.at(xk))\n",
        "    if k < N-1:\n",
        "      controls.append(result.at(uk))\n",
        "\n",
        "\n",
        "states = np.array(states)\n",
        "controls = np.array(controls)\n",
        "\n",
        "# Plotting\n",
        "plt.figure(figsize=(10, 4))\n",
        "plt.plot(states[:, 0], label='Position')\n",
        "plt.plot(states[:, 1], label='Velocity')\n",
        "plt.title('GTSAM-solved State Trajectories')\n",
        "plt.legend()\n",
        "states\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HyEtrvX1bigY"
      },
      "source": [
        "###  Visualizing the Factor Graph and Variable Elimination (Bayes Net)\n",
        "\n",
        "####  Initial Factor Graph\n",
        "\n",
        "The first plot generated by `graphviz.Source(graph.dot())` visualizes the **factor graph** representation of the LQR problem.\n",
        "\n",
        "- **Ellipses** represent variable nodes: states $ x_k $ and controls $ u_k $.\n",
        "- **Black dots** are factor nodes connecting variables — each represents either:\n",
        "  - A cost factor (state or control penalty),\n",
        "  - A dynamics constraint,\n",
        "  - Or a prior (for $ x_0 $).\n",
        "\n",
        "The structure reflects the temporal nature of the problem:\n",
        "- Each state $ x_k $ is connected to its predecessor and successor via dynamics,\n",
        "- Each control $ u_k $ connects two consecutive states,\n",
        "- Cost factors attach individually to each $ x_k $ and $ u_k $.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 102
        },
        "id": "ggr3ukXoMw06",
        "outputId": "7635d384-d4e3-4565-9639-541f9bba23de"
      },
      "outputs": [
        {
          "data": {
            "image/svg+xml": [
              "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
              "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
              " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
              "<!-- Generated by graphviz version 2.43.0 (0)\n",
              " -->\n",
              "<!-- Title: %3 Pages: 1 -->\n",
              "<svg width=\"360pt\" height=\"61pt\"\n",
              " viewBox=\"0.00 0.00 360.00 60.92\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
              "<g id=\"graph0\" class=\"graph\" transform=\"scale(0.73 0.73) rotate(0) translate(4 79.6)\">\n",
              "<title>%3</title>\n",
              "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-79.6 490,-79.6 490,4 -4,4\"/>\n",
              "<!-- var8430738502437568512 -->\n",
              "<g id=\"node1\" class=\"node\">\n",
              "<title>var8430738502437568512</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"27\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">u0</text>\n",
              "</g>\n",
              "<!-- factor2 -->\n",
              "<g id=\"node10\" class=\"node\">\n",
              "<title>factor2</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"27\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&#45;factor2 -->\n",
              "<g id=\"edge3\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&#45;factor2</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M27,-39.58C27,-25.79 27,-7.97 27,-3.73\"/>\n",
              "</g>\n",
              "<!-- factor3 -->\n",
              "<g id=\"node11\" class=\"node\">\n",
              "<title>factor3</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"77\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&#45;factor3 -->\n",
              "<g id=\"edge6\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&#45;factor3</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M40.66,-41.9C53.62,-27.96 71.72,-8.48 76.05,-3.83\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513 -->\n",
              "<g id=\"node2\" class=\"node\">\n",
              "<title>var8430738502437568513</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"243\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">u1</text>\n",
              "</g>\n",
              "<!-- factor5 -->\n",
              "<g id=\"node13\" class=\"node\">\n",
              "<title>factor5</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"254\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&#45;factor5 -->\n",
              "<g id=\"edge8\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&#45;factor5</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M246.48,-39.58C249.3,-25.79 252.94,-7.97 253.81,-3.73\"/>\n",
              "</g>\n",
              "<!-- factor6 -->\n",
              "<g id=\"node14\" class=\"node\">\n",
              "<title>factor6</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"232\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&#45;factor6 -->\n",
              "<g id=\"edge11\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&#45;factor6</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M239.52,-39.58C236.7,-25.79 233.06,-7.97 232.19,-3.73\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514 -->\n",
              "<g id=\"node3\" class=\"node\">\n",
              "<title>var8430738502437568514</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"387\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"387\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">u2</text>\n",
              "</g>\n",
              "<!-- factor8 -->\n",
              "<g id=\"node16\" class=\"node\">\n",
              "<title>factor8</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"398\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&#45;factor8 -->\n",
              "<g id=\"edge13\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&#45;factor8</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M390.48,-39.58C393.3,-25.79 396.94,-7.97 397.81,-3.73\"/>\n",
              "</g>\n",
              "<!-- factor9 -->\n",
              "<g id=\"node17\" class=\"node\">\n",
              "<title>factor9</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"376\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&#45;factor9 -->\n",
              "<g id=\"edge16\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&#45;factor9</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M383.52,-39.58C380.7,-25.79 377.06,-7.97 376.19,-3.73\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320 -->\n",
              "<g id=\"node4\" class=\"node\">\n",
              "<title>var8646911284551352320</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"99\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
              "</g>\n",
              "<!-- factor0 -->\n",
              "<g id=\"node8\" class=\"node\">\n",
              "<title>factor0</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"99\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320&#45;&#45;factor0 -->\n",
              "<g id=\"edge1\" class=\"edge\">\n",
              "<title>var8646911284551352320&#45;&#45;factor0</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M99,-39.58C99,-25.79 99,-7.97 99,-3.73\"/>\n",
              "</g>\n",
              "<!-- factor1 -->\n",
              "<g id=\"node9\" class=\"node\">\n",
              "<title>factor1</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"121\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320&#45;&#45;factor1 -->\n",
              "<g id=\"edge2\" class=\"edge\">\n",
              "<title>var8646911284551352320&#45;&#45;factor1</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M105.72,-40.17C111.38,-26.32 118.83,-8.1 120.61,-3.76\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320&#45;&#45;factor3 -->\n",
              "<g id=\"edge5\" class=\"edge\">\n",
              "<title>var8646911284551352320&#45;&#45;factor3</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M92.28,-40.17C86.62,-26.32 79.17,-8.1 77.39,-3.76\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321 -->\n",
              "<g id=\"node5\" class=\"node\">\n",
              "<title>var8646911284551352321</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"171\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
              "</g>\n",
              "<!-- var8646911284551352321&#45;&#45;factor3 -->\n",
              "<g id=\"edge4\" class=\"edge\">\n",
              "<title>var8646911284551352321&#45;&#45;factor3</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M151.07,-45.19C126.06,-30.88 85.55,-7.69 78.17,-3.47\"/>\n",
              "</g>\n",
              "<!-- factor4 -->\n",
              "<g id=\"node12\" class=\"node\">\n",
              "<title>factor4</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"171\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321&#45;&#45;factor4 -->\n",
              "<g id=\"edge7\" class=\"edge\">\n",
              "<title>var8646911284551352321&#45;&#45;factor4</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M171,-39.58C171,-25.79 171,-7.97 171,-3.73\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321&#45;&#45;factor6 -->\n",
              "<g id=\"edge10\" class=\"edge\">\n",
              "<title>var8646911284551352321&#45;&#45;factor6</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M186.71,-42.75C202.52,-28.8 225.34,-8.68 230.8,-3.86\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322 -->\n",
              "<g id=\"node6\" class=\"node\">\n",
              "<title>var8646911284551352322</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"315\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"315\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
              "</g>\n",
              "<!-- var8646911284551352322&#45;&#45;factor6 -->\n",
              "<g id=\"edge9\" class=\"edge\">\n",
              "<title>var8646911284551352322&#45;&#45;factor6</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M296.17,-44.39C273.97,-30.01 239.31,-7.54 233,-3.45\"/>\n",
              "</g>\n",
              "<!-- factor7 -->\n",
              "<g id=\"node15\" class=\"node\">\n",
              "<title>factor7</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"315\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322&#45;&#45;factor7 -->\n",
              "<g id=\"edge12\" class=\"edge\">\n",
              "<title>var8646911284551352322&#45;&#45;factor7</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M315,-39.58C315,-25.79 315,-7.97 315,-3.73\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322&#45;&#45;factor9 -->\n",
              "<g id=\"edge15\" class=\"edge\">\n",
              "<title>var8646911284551352322&#45;&#45;factor9</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M330.71,-42.75C346.52,-28.8 369.34,-8.68 374.8,-3.86\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352323 -->\n",
              "<g id=\"node7\" class=\"node\">\n",
              "<title>var8646911284551352323</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"459\" cy=\"-57.6\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"459\" y=\"-53.9\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
              "</g>\n",
              "<!-- var8646911284551352323&#45;&#45;factor9 -->\n",
              "<g id=\"edge14\" class=\"edge\">\n",
              "<title>var8646911284551352323&#45;&#45;factor9</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M440.17,-44.39C417.97,-30.01 383.31,-7.54 377,-3.45\"/>\n",
              "</g>\n",
              "<!-- factor10 -->\n",
              "<g id=\"node18\" class=\"node\">\n",
              "<title>factor10</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"448\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352323&#45;&#45;factor10 -->\n",
              "<g id=\"edge17\" class=\"edge\">\n",
              "<title>var8646911284551352323&#45;&#45;factor10</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M455.52,-39.58C452.7,-25.79 449.06,-7.97 448.19,-3.73\"/>\n",
              "</g>\n",
              "<!-- factor11 -->\n",
              "<g id=\"node19\" class=\"node\">\n",
              "<title>factor11</title>\n",
              "<ellipse fill=\"black\" stroke=\"black\" cx=\"470\" cy=\"-1.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352323&#45;&#45;factor11 -->\n",
              "<g id=\"edge18\" class=\"edge\">\n",
              "<title>var8646911284551352323&#45;&#45;factor11</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M462.48,-39.58C465.3,-25.79 468.94,-7.97 469.81,-3.73\"/>\n",
              "</g>\n",
              "</g>\n",
              "</svg>\n"
            ],
            "text/plain": [
              "<graphviz.sources.Source at 0x7e2318659510>"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "graphviz.Source(graph.dot())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Re3c12XjdCfu"
      },
      "source": [
        "###  From Factor Graph to Bayes Net via Variable Elimination\n",
        "\n",
        "After constructing a Gaussian factor graph to represent the finite-horizon LQR problem, we can use GTSAM's `eliminateSequential(ordering)` to transform the factor graph into a **Bayes Net** — a directed acyclic graphical model representing the posterior distribution over all variables in a factored form.\n",
        "\n",
        "The following plots will visualize the **Bayes Net** obtained using the elimination order:\n",
        "\n",
        "$$\n",
        "\\text{Ordering} = [x_0, u_0, x_1, u_1, x_2, u_2, x_3]\n",
        "$$\n",
        "\n",
        "---\n",
        "\n",
        "###  Structure of the Bayes Net\n",
        "\n",
        "Each node in the plot represents a **variable**, and each arrow represents a **dependency** in a conditional Gaussian distribution arising from the variable elimination process. For example:\n",
        "\n",
        "- $ x_2 \\leftarrow \\{x_3, u_2\\} $  \n",
        "  ⇒ $ p(x_2 \\mid x_3, u_2) $\n",
        "- $ u_1 \\leftarrow \\{x_2\\} $  \n",
        "  ⇒ $ p(u_1 \\mid x_2) $\n",
        "  \n",
        "This structure means that once we solve for the \"latest\" variable (e.g., $ x_3 $), we can **back-substitute** through the DAG to recover the entire trajectory.\n",
        "\n",
        "---\n",
        "\n",
        "### What Does Elimination Do?\n",
        "\n",
        "In mathematical terms, elimination rewrites the joint Gaussian density:\n",
        "\n",
        "$\n",
        "p(x_0, u_0, x_1, u_1, \\dots, x_N) \\propto \\prod_i \\phi_i(x_i, u_i)\n",
        "$\n",
        "\n",
        "as a sequence of **conditionals**:\n",
        "\n",
        "$\n",
        "p(x_0, u_0, x_1, \\dots) = p(x_3) \\cdot p(x_2 \\mid x_3, u_2) \\cdot \\dots \\cdot p(x_0)\n",
        "$\n",
        "\n",
        "Each conditional is Gaussian:\n",
        "\n",
        "$\n",
        "x_k = A_k \\cdot \\text{Parents}(x_k) + \\eta_k, \\quad \\eta_k \\sim \\mathcal{N}(0, \\Sigma_k)\n",
        "$\n",
        "\n",
        "This allows us to represent the problem as a **triangular system**, suitable for efficient solving.\n",
        "\n",
        "---\n",
        "\n",
        "###  Effect of Elimination Ordering\n",
        "\n",
        "####  Will the result (optimal trajectory) change if we change the elimination ordering?\n",
        "\n",
        "No — the **MAP solution** is invariant under ordering for **linear-Gaussian systems**, because the solution to:\n",
        "\n",
        "$$\n",
        "\\min_x \\|Ax - b\\|^2\n",
        "$$\n",
        "\n",
        "is unique and does **not** depend on the factorization path.\n",
        "\n",
        "####  What changes?\n",
        "\n",
        "1. **Structure of the Bayes Net**:  \n",
        "   The set of conditionals and their dependencies change. E.g., if we eliminate all controls first, each $ u_k $ may depend on more than just $ x_k $, introducing wider conditionals.\n",
        "\n",
        "2. **Computational efficiency**:  \n",
        "   Different orderings lead to different **fill-in patterns** during matrix factorization (e.g., Cholesky). A poor ordering increases memory and time.\n",
        "\n",
        "3. **Sparsity pattern**:  \n",
        "   An optimized ordering (e.g., using COLAMD or METIS) maintains sparse Jacobians and minimal cross-dependencies.\n",
        "\n",
        "4. **Numerical stability** (nonlinear problems):  \n",
        "   Poor ordering may lead to ill-conditioned systems or less stable Gauss-Newton updates.\n",
        "\n",
        "---\n",
        "\n",
        "###  Analogy to QR/Cholesky\n",
        "\n",
        "The process is equivalent to **structured Gaussian elimination** or **QR factorization** of the system $ Ax = b $, where variable elimination corresponds to pivoting rows and columns of the matrix in a specific order.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FZ1zz1ozM0RP",
        "outputId": "429d71de-5c73-47ba-e9af-9ad2606af0c4"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Ordering: Position 0: x3, u2, x2, u1, x1, u0, x0\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ordering = gtsam.Ordering()\n",
        "for k in range(N):\n",
        "  ordering.push_back(X(k))\n",
        "  if k < N-1:\n",
        "    ordering.push_back(U(k))\n",
        "\n",
        "\n",
        "# Generate a reversed ordering from above order\n",
        "ordering_reversed = gtsam.Ordering()\n",
        "for k in range(N-1, -1, -1):\n",
        "  ordering_reversed.push_back(X(k))\n",
        "  if k > 0:\n",
        "    ordering_reversed.push_back(U(k-1))\n",
        "\n",
        "\n",
        "\n",
        "print(\"Ordering:\", ordering_reversed)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 175
        },
        "id": "8ZT2qtwyM7o_",
        "outputId": "63a0a997-0196-4572-b80a-b62af971d88c"
      },
      "outputs": [
        {
          "data": {
            "image/svg+xml": [
              "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
              "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
              " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
              "<!-- Generated by graphviz version 2.43.0 (0)\n",
              " -->\n",
              "<!-- Title: %3 Pages: 1 -->\n",
              "<svg width=\"278pt\" height=\"116pt\"\n",
              " viewBox=\"0.00 0.00 278.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
              "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
              "<title>%3</title>\n",
              "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 274,-112 274,4 -4,4\"/>\n",
              "<!-- var8430738502437568512 -->\n",
              "<g id=\"node1\" class=\"node\">\n",
              "<title>var8430738502437568512</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u0</text>\n",
              "</g>\n",
              "<!-- var8646911284551352321 -->\n",
              "<g id=\"node5\" class=\"node\">\n",
              "<title>var8646911284551352321</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&gt;var8646911284551352321 -->\n",
              "<g id=\"edge2\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&gt;var8646911284551352321</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M42.27,-74.73C52.2,-64.8 65.32,-51.68 76.44,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"79.16,-42.79 83.75,-33.25 74.21,-37.84 79.16,-42.79\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513 -->\n",
              "<g id=\"node2\" class=\"node\">\n",
              "<title>var8430738502437568513</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u1</text>\n",
              "</g>\n",
              "<!-- var8646911284551352322 -->\n",
              "<g id=\"node6\" class=\"node\">\n",
              "<title>var8646911284551352322</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"171\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&gt;var8646911284551352322 -->\n",
              "<g id=\"edge5\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&gt;var8646911284551352322</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M114.27,-74.73C124.2,-64.8 137.32,-51.68 148.44,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"151.16,-42.79 155.75,-33.25 146.21,-37.84 151.16,-42.79\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514 -->\n",
              "<g id=\"node3\" class=\"node\">\n",
              "<title>var8430738502437568514</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"171\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u2</text>\n",
              "</g>\n",
              "<!-- var8646911284551352323 -->\n",
              "<g id=\"node7\" class=\"node\">\n",
              "<title>var8646911284551352323</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"243\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&gt;var8646911284551352323 -->\n",
              "<g id=\"edge8\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&gt;var8646911284551352323</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M186.27,-74.73C196.2,-64.8 209.32,-51.68 220.44,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"223.16,-42.79 227.75,-33.25 218.21,-37.84 223.16,-42.79\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320 -->\n",
              "<g id=\"node4\" class=\"node\">\n",
              "<title>var8646911284551352320</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
              "</g>\n",
              "<!-- var8646911284551352320&#45;&gt;var8430738502437568512 -->\n",
              "<g id=\"edge1\" class=\"edge\">\n",
              "<title>var8646911284551352320&#45;&gt;var8430738502437568512</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M27,-36.17C27,-43.87 27,-53.03 27,-61.58\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"23.5,-61.59 27,-71.59 30.5,-61.59 23.5,-61.59\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320&#45;&gt;var8646911284551352321 -->\n",
              "<g id=\"edge3\" class=\"edge\">\n",
              "<title>var8646911284551352320&#45;&gt;var8646911284551352321</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M54.22,-18C56.64,-18 59.11,-18 61.6,-18\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"61.74,-21.5 71.74,-18 61.74,-14.5 61.74,-21.5\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321&#45;&gt;var8430738502437568513 -->\n",
              "<g id=\"edge4\" class=\"edge\">\n",
              "<title>var8646911284551352321&#45;&gt;var8430738502437568513</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M99,-36.17C99,-43.87 99,-53.03 99,-61.58\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"95.5,-61.59 99,-71.59 102.5,-61.59 95.5,-61.59\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321&#45;&gt;var8646911284551352322 -->\n",
              "<g id=\"edge6\" class=\"edge\">\n",
              "<title>var8646911284551352321&#45;&gt;var8646911284551352322</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M126.22,-18C128.64,-18 131.11,-18 133.6,-18\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"133.74,-21.5 143.74,-18 133.74,-14.5 133.74,-21.5\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322&#45;&gt;var8430738502437568514 -->\n",
              "<g id=\"edge7\" class=\"edge\">\n",
              "<title>var8646911284551352322&#45;&gt;var8430738502437568514</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M171,-36.17C171,-43.87 171,-53.03 171,-61.58\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"167.5,-61.59 171,-71.59 174.5,-61.59 167.5,-61.59\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322&#45;&gt;var8646911284551352323 -->\n",
              "<g id=\"edge9\" class=\"edge\">\n",
              "<title>var8646911284551352322&#45;&gt;var8646911284551352323</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M198.22,-18C200.64,-18 203.11,-18 205.6,-18\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"205.74,-21.5 215.74,-18 205.74,-14.5 205.74,-21.5\"/>\n",
              "</g>\n",
              "</g>\n",
              "</svg>\n"
            ],
            "text/plain": [
              "<gtbook.display.show at 0x7e22fb05f3d0>"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "dag = graph.eliminateSequential(ordering_reversed)\n",
        "show(dag, hints={\"u\":1, \"x\":0})"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 448
        },
        "id": "0MyYkdtkYdus",
        "outputId": "593a12b3-bfa7-4b4b-fd63-2e0091eda231"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7e23182eb450>"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAGdCAYAAADJ82znAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIBlJREFUeJzt3XtwVPX9//HXJsDGMslGS0iyEm4qoFyCosQgFikpIWWQ0FYxQ8tFwA7f0JFfqtU4SkCdxnrDKvmC7RiCQ5XLDAZHmdgQBYrc5JKRUMuQNJAwZIMwJktiDTR7fn/0y9qVbHTJ+SRLfD5mPjOccz6fT97nxH159iT7icOyLEsAYEBEVxcAoPsiYAAYQ8AAMIaAAWAMAQPAGAIGgDEEDABjCBgAxvTo6gLs4PP5dPr0aUVHR8vhcHR1OUC3YFmWzp8/L7fbrYiIK7sX6RYBc/r0aSUlJXV1GUC3VFtbq379+l3R2G4RMNHR0f/5xy9HS70iOzxfw8t5HZ4DuNp5vV+qf/8Hvn59XYFuETD+t0W9IqVeHT+lmJjeHZ4D6C468tiBh7wAjCFgABhjLGAKCgo0cOBARUVFKSUlRfv372+3/6ZNmzRs2DBFRUVp5MiR2rp1q6nSAHQSIwGzYcMG5eTkKC8vT4cOHVJycrLS09N15syZNvvv3r1bWVlZmj9/vg4fPqzMzExlZmaqoqLCRHkAOonDxIJTKSkpuuOOO7Ry5UpJ//k9laSkJP3mN7/R448/fln/mTNnqrm5We+9955/35133qnRo0dr9erV3/r1vF6vXC6X9OAYWx7y+v739x2eA7jaeb3Nio29V42NjYqJibmiOWy/g7lw4YIOHjyotLS0r79IRITS0tK0Z8+eNsfs2bMnoL8kpaenB+3f0tIir9cb0ACEH9sD5uzZs2ptbVV8fHzA/vj4eHk8njbHeDyekPrn5+fL5XL5G79kB4Snq/KnSLm5uWpsbPS32trari4JQBts/0W7Pn36KDIyUvX19QH76+vrlZCQ0OaYhISEkPo7nU45nU57CgZgjO13ML169dKYMWNUVlbm3+fz+VRWVqbU1NQ2x6Smpgb0l6TS0tKg/QFcHYx8VCAnJ0dz5szR7bffrrFjx+qVV15Rc3Oz5s2bJ0maPXu2rr/+euXn50uSHn74YU2YMEEvvfSSpk6dqvXr1+vAgQP605/+ZKI8AJ3ESMDMnDlTn3/+uZYuXSqPx6PRo0erpKTE/yC3pqYm4OPf48aN01tvvaUnn3xSTzzxhG666SYVFxdrxIgRJsoD0EmM/B5MZ+P3YAD7heXvwQDAJQQMAGO6xXowlzS8nGfLWi4R//OEDdV8jbdc+L7iDgaAMQQMAGMIGADGEDAAjCFgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGdKs1ee1i9xq6dq7xy/q+uJpwBwPAGAIGgDEEDABjCBgAxhAwAIwhYAAYQ8AAMIaAAWAMAQPAGAIGgDEEDABjCBgAxhAwAIwhYAAYQ8AAMMb2gMnPz9cdd9yh6Oho9e3bV5mZmTp27Fi7Y4qKiuRwOAJaVFSU3aUB6GS2B8yOHTuUnZ2tvXv3qrS0VBcvXtTkyZPV3Nzc7riYmBjV1dX528mTJ+0uDUAns31Fu5KSkoDtoqIi9e3bVwcPHtSPfvSjoOMcDocSEhLsLgdAFzK+ZGZjY6Mk6brrrmu3X1NTkwYMGCCfz6fbbrtNv//97zV8+PA2+7a0tKilpcW/7fV6JUmbYv+ffqDIDtec9UVeh+f4b3Yuc2nn8psSS3DCLKMPeX0+n5YsWaK77rpLI0aMCNpv6NChKiws1JYtW7Ru3Tr5fD6NGzdOp06darN/fn6+XC6XvyUlJZk6BQAdYDRgsrOzVVFRofXr17fbLzU1VbNnz9bo0aM1YcIEbd68WXFxcXr99dfb7J+bm6vGxkZ/q62tNVE+gA4y9hZp8eLFeu+997Rz507169cvpLE9e/bUrbfeqsrKyjaPO51OOZ1OO8oEYJDtdzCWZWnx4sV655139OGHH2rQoEEhz9Ha2qojR44oMTHR7vIAdCLb72Cys7P11ltvacuWLYqOjpbH45EkuVwuXXPNNZKk2bNn6/rrr1d+fr4k6emnn9add96pG2+8UQ0NDXrhhRd08uRJLViwwO7yAHQi2wNm1apVkqR77rknYP+aNWs0d+5cSVJNTY0iIr6+efriiy+0cOFCeTweXXvttRozZox2796tW265xe7yAHQi2wPGsqxv7bN9+/aA7RUrVmjFihV2lwKgi/FZJADGEDAAjCFgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwxviavJ3pvoYVionp3eF53o5YZEM1Zvh8q2ydz841flnfF9/EHQwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMZ0qyUz7ZJl87KUdrJziUuJZS5hFncwAIwhYAAYQ8AAMIaAAWAMAQPAGNsDZtmyZXI4HAFt2LBh7Y7ZtGmThg0bpqioKI0cOVJbt261uywAXcDIHczw4cNVV1fnb7t27Qrad/fu3crKytL8+fN1+PBhZWZmKjMzUxUVFSZKA9CJjARMjx49lJCQ4G99+vQJ2vePf/yjpkyZokcffVQ333yznnnmGd12221auXKlidIAdCIjAXP8+HG53W4NHjxYs2bNUk1NTdC+e/bsUVpaWsC+9PR07dmzx0RpADqR7b/Jm5KSoqKiIg0dOlR1dXVavny57r77blVUVCg6Ovqy/h6PR/Hx8QH74uPj5fF4gn6NlpYWtbS0+Le9Xq99JwDANrYHTEZGhv/fo0aNUkpKigYMGKCNGzdq/vz5tnyN/Px8LV++3Ja5AJhj/MfUsbGxGjJkiCorK9s8npCQoPr6+oB99fX1SkhICDpnbm6uGhsb/a22ttbWmgHYw3jANDU1qaqqSomJiW0eT01NVVlZWcC+0tJSpaamBp3T6XQqJiYmoAEIP7YHzCOPPKIdO3boxIkT2r17t2bMmKHIyEhlZWVJkmbPnq3c3Fx//4cfflglJSV66aWX9I9//EPLli3TgQMHtHjxYrtLA9DJbH8Gc+rUKWVlZencuXOKi4vT+PHjtXfvXsXFxUmSampqFBHxda6NGzdOb731lp588kk98cQTuummm1RcXKwRI0bYXRqATuawLMvq6iI6yuv1yuVyqaHhXcXE9O7qcoxiPRh0Fq+3WbGx96qxsfGKH0PwWSQAxhAwAIwhYAAYw5q8Vxm7n5nY+UyH5zn4Ju5gABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwhiUzv+fsXOaSP6mCb+IOBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgjO0BM3DgQDkcjstadnZ2m/2Lioou6xsVFWV3WQC6gO0LTn3yySdqbW31b1dUVOgnP/mJ7rvvvqBjYmJidOzYMf+2w+GwuywAXcD2gImLiwvYfu6553TDDTdowoQJQcc4HA4lJCTYXQqALmb0GcyFCxe0bt06Pfjgg+3elTQ1NWnAgAFKSkrS9OnTdfTo0XbnbWlpkdfrDWgAwo/RNXmLi4vV0NCguXPnBu0zdOhQFRYWatSoUWpsbNSLL76ocePG6ejRo+rXr1+bY/Lz87V8+XJDVYe3tyMW2Tpf1hd5ts1l9xq6dq7xy/q+XcPoHcwbb7yhjIwMud3uoH1SU1M1e/ZsjR49WhMmTNDmzZsVFxen119/PeiY3NxcNTY2+lttba2J8gF0kLE7mJMnT2rbtm3avHlzSON69uypW2+9VZWVlUH7OJ1OOZ3OjpYIwDBjdzBr1qxR3759NXXq1JDGtba26siRI0pMTDRUGYDOYiRgfD6f1qxZozlz5qhHj8CbpNmzZys3N9e//fTTT+uvf/2r/vnPf+rQoUP65S9/qZMnT2rBggUmSgPQiYy8Rdq2bZtqamr04IMPXnaspqZGERFf59oXX3yhhQsXyuPx6Nprr9WYMWO0e/du3XLLLSZKA9CJjATM5MmTZVlWm8e2b98esL1ixQqtWLHCRBkAuhifRQJgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAY4wumQn7ZflW2Tqf3Utw2sln47naufymxBKc3xV3MACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGNXm/5+xe49dOdq6jyxq6XYM7GADGEDAAjCFgABhDwAAwhoABYEzIAbNz505NmzZNbrdbDodDxcXFAccty9LSpUuVmJioa665RmlpaTp+/Pi3zltQUKCBAwcqKipKKSkp2r9/f6ilAQgzIQdMc3OzkpOTVVBQ0Obx559/Xq+++qpWr16tffv2qXfv3kpPT9dXX30VdM4NGzYoJydHeXl5OnTokJKTk5Wenq4zZ86EWh6AMBJywGRkZOjZZ5/VjBkzLjtmWZZeeeUVPfnkk5o+fbpGjRqlN998U6dPn77sTue/vfzyy1q4cKHmzZunW265RatXr9YPfvADFRYWhloegDBi6zOY6upqeTwepaWl+fe5XC6lpKRoz549bY65cOGCDh48GDAmIiJCaWlpQce0tLTI6/UGNADhx9aA8Xg8kqT4+PiA/fHx8f5j33T27Fm1traGNCY/P18ul8vfkpKSbKgegN2uyp8i5ebmqrGx0d9qa2u7uiQAbbA1YBISEiRJ9fX1Afvr6+v9x76pT58+ioyMDGmM0+lUTExMQAMQfmwNmEGDBikhIUFlZWX+fV6vV/v27VNqamqbY3r16qUxY8YEjPH5fCorKws6BsDVIeRPUzc1NamystK/XV1drfLycl133XXq37+/lixZomeffVY33XSTBg0apKeeekput1uZmZn+MZMmTdKMGTO0ePFiSVJOTo7mzJmj22+/XWPHjtUrr7yi5uZmzZs3r+NnCKDLhBwwBw4c0MSJE/3bOTk5kqQ5c+aoqKhIv/vd79Tc3KyHHnpIDQ0NGj9+vEpKShQVFeUfU1VVpbNnz/q3Z86cqc8//1xLly6Vx+PR6NGjVVJSctmDXwBXF4dlWVZXF9FRXq9XLpdLDQ3vKiamd1eXA5uwHkzX8nqbFRt7rxobG6/4OedV+VMkAFcHAgaAMSyZibBl59saO99uSbzl+q64gwFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxrMmL7wW719DlT6p8N9zBADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMaEHDA7d+7UtGnT5Ha75XA4VFxc7D928eJFPfbYYxo5cqR69+4tt9ut2bNn6/Tp0+3OuWzZMjkcjoA2bNiwkE8GQHgJOWCam5uVnJysgoKCy459+eWXOnTokJ566ikdOnRImzdv1rFjx3Tvvfd+67zDhw9XXV2dv+3atSvU0gCEmZAXnMrIyFBGRkabx1wul0pLSwP2rVy5UmPHjlVNTY369+8fvJAePZSQkBBqOQDCmPFnMI2NjXI4HIqNjW233/Hjx+V2uzV48GDNmjVLNTU1pksDYJjRJTO/+uorPfbYY8rKylJMTEzQfikpKSoqKtLQoUNVV1en5cuX6+6771ZFRYWio6Mv69/S0qKWlhb/ttfrNVI/utbbEYtsmyvrizzb5pLsXebSzuU3pfBagtNYwFy8eFH333+/LMvSqlWr2u3732+5Ro0apZSUFA0YMEAbN27U/PnzL+ufn5+v5cuX214zAHsZeYt0KVxOnjyp0tLSdu9e2hIbG6shQ4aosrKyzeO5ublqbGz0t9raWjvKBmAz2wPmUrgcP35c27Zt0w9/+MOQ52hqalJVVZUSExPbPO50OhUTExPQAISfkAOmqalJ5eXlKi8vlyRVV1ervLxcNTU1unjxon7xi1/owIED+stf/qLW1lZ5PB55PB5duHDBP8ekSZO0cuVK//YjjzyiHTt26MSJE9q9e7dmzJihyMhIZWVldfwMAXSZkJ/BHDhwQBMnTvRv5+TkSJLmzJmjZcuW6d1335UkjR49OmDcRx99pHvuuUeSVFVVpbNnz/qPnTp1SllZWTp37pzi4uI0fvx47d27V3FxcaGWByCMhBww99xzjyzLCnq8vWOXnDhxImB7/fr1oZYB4CrAZ5EAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxRtfkBToiy9f+UquhsHN9X7v5bDxPycY1fi/8u8NTcAcDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxLJmJ7wU7l9+0m21LXP4f3//+3pZ5vN5mxRbe26E5uIMBYAwBA8AYAgaAMQQMAGMIGADGhBwwO3fu1LRp0+R2u+VwOFRcXBxwfO7cuXI4HAFtypQp3zpvQUGBBg4cqKioKKWkpGj//v2hlgYgzIQcMM3NzUpOTlZBQUHQPlOmTFFdXZ2/vf322+3OuWHDBuXk5CgvL0+HDh1ScnKy0tPTdebMmVDLAxBGQv49mIyMDGVkZLTbx+l0KiEh4TvP+fLLL2vhwoWaN2+eJGn16tV6//33VVhYqMcffzzUEgGECSPPYLZv366+fftq6NChWrRokc6dOxe074ULF3Tw4EGlpaV9XVREhNLS0rRnz542x7S0tMjr9QY0AOHH9oCZMmWK3nzzTZWVlekPf/iDduzYoYyMDLW2trbZ/+zZs2ptbVV8fHzA/vj4eHk8njbH5Ofny+Vy+VtSUpLdpwHABrZ/VOCBBx7w/3vkyJEaNWqUbrjhBm3fvl2TJk2y5Wvk5uYqJyfHv+31egkZIAwZ/zH14MGD1adPH1VWVrZ5vE+fPoqMjFR9fX3A/vr6+qDPcZxOp2JiYgIagPBjPGBOnTqlc+fOKTExsc3jvXr10pgxY1RWVubf5/P5VFZWptTUVNPlATAo5IBpampSeXm5ysvLJUnV1dUqLy9XTU2Nmpqa9Oijj2rv3r06ceKEysrKNH36dN14441KT0/3zzFp0iStXLnSv52Tk6M///nPWrt2rT777DMtWrRIzc3N/p8qAbg6hfwM5sCBA5o4caJ/+9KzkDlz5mjVqlX69NNPtXbtWjU0NMjtdmvy5Ml65pln5HQ6/WOqqqp09uxZ//bMmTP1+eefa+nSpfJ4PBo9erRKSkoue/AL4OrisCzL6uoiOsrr9crlcqmh4V3FxPTu6nKAkIT1ejCx96qxsfGKn3PyWSQAxhAwAIwhYAAYw5q8QBez65nJJbY907nw7w5PwR0MAGMIGADGEDAAjCFgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGsGQm0MXC+s+WFN7boTm4gwFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGEPAADCGgAFgDAEDwBgCBoAxBAwAYwgYAMYQMACMIWAAGBNywOzcuVPTpk2T2+2Ww+FQcXFxwHGHw9Fme+GFF4LOuWzZssv6Dxs2LOSTARBeQg6Y5uZmJScnq6CgoM3jdXV1Aa2wsFAOh0M///nP2513+PDhAeN27doVamkAwkzIK9plZGQoIyMj6PGEhISA7S1btmjixIkaPHhw+4X06HHZWABXN6NLZtbX1+v999/X2rVrv7Xv8ePH5Xa7FRUVpdTUVOXn56t///5t9m1paVFLS4t/2+v12lYz0NnsWuLyEtuW4Lzw7w5PYfQh79q1axUdHa2f/exn7fZLSUlRUVGRSkpKtGrVKlVXV+vuu+/W+fPn2+yfn58vl8vlb0lJSSbKB9BBRgOmsLBQs2bNUlRUVLv9MjIydN9992nUqFFKT0/X1q1b1dDQoI0bN7bZPzc3V42Njf5WW1tronwAHWTsLdLf/vY3HTt2TBs2bAh5bGxsrIYMGaLKyso2jzudTjmdzo6WCMAwY3cwb7zxhsaMGaPk5OSQxzY1NamqqkqJiYkGKgPQWUIOmKamJpWXl6u8vFySVF1drfLyctXU1Pj7eL1ebdq0SQsWLGhzjkmTJmnlypX+7UceeUQ7duzQiRMntHv3bs2YMUORkZHKysoKtTwAYSTkt0gHDhzQxIkT/ds5OTmSpDlz5qioqEiStH79elmWFTQgqqqqdPbsWf/2qVOnlJWVpXPnzikuLk7jx4/X3r17FRcXF2p5AMKIw7Isq6uL6Civ1yuXy6WGhncVE9O7q8sBupStP6YuPKjGxkbFxMRcWS32VAIAlyNgABhDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwxuiSmZ3l0sepvN4vu7gSIAzYsNTlf+ZplfT16+tKdIuAubS0Zv/+D3RxJUD3c/78eblcrisa2y0+Te3z+XT69GlFR0fL4XAE7ef1epWUlKTa2tor/nRoV7vaz+Fqr1/6/pyDZVk6f/683G63IiKu7GlKt7iDiYiIUL9+/b5z/5iYmKv2P4xLrvZzuNrrl74f53Cldy6X8JAXgDEEDABjvlcB43Q6lZeXd1X/RYKr/Ryu9volziEU3eIhL4Dw9L26gwHQuQgYAMYQMACMIWAAGNPtAqagoEADBw5UVFSUUlJStH///nb7b9q0ScOGDVNUVJRGjhyprVu3dlKll8vPz9cdd9yh6Oho9e3bV5mZmTp27Fi7Y4qKiuRwOAJaVFRUJ1V8uWXLll1Wz7Bhw9odE07fg4EDB15Wv8PhUHZ2dpv9w+H679y5U9OmTZPb7ZbD4VBxcXHAccuytHTpUiUmJuqaa65RWlqajh8//q3zhvpaaku3CpgNGzYoJydHeXl5OnTokJKTk5Wenq4zZ8602X/37t3KysrS/PnzdfjwYWVmZiozM1MVFRWdXPl/7NixQ9nZ2dq7d69KS0t18eJFTZ48Wc3Nze2Oi4mJUV1dnb+dPHmykypu2/DhwwPq2bVrV9C+4fY9+OSTTwJqLy0tlSTdd999Qcd09fVvbm5WcnKyCgoK2jz+/PPP69VXX9Xq1au1b98+9e7dW+np6frqq6+CzhnqaykoqxsZO3aslZ2d7d9ubW213G63lZ+f32b/+++/35o6dWrAvpSUFOvXv/610Tq/qzNnzliSrB07dgTts2bNGsvlcnVeUd8iLy/PSk5O/s79w/178PDDD1s33HCD5fP52jwebtdfkvXOO+/4t30+n5WQkGC98MIL/n0NDQ2W0+m03n777aDzhPpaCqbb3MFcuHBBBw8eVFpamn9fRESE0tLStGfPnjbH7NmzJ6C/JKWnpwft39kaGxslSdddd127/ZqamjRgwAAlJSVp+vTpOnr0aGeUF9Tx48fldrs1ePBgzZo1SzU1NUH7hvP34MKFC1q3bp0efPDBdj9EG27X/79VV1fL4/EEXGOXy6WUlJSg1/hKXkvBdJuAOXv2rFpbWxUfHx+wPz4+Xh6Pp80xHo8npP6dyefzacmSJbrrrrs0YsSIoP2GDh2qwsJCbdmyRevWrZPP59O4ceN06tSpTqz2aykpKSoqKlJJSYlWrVql6upq3X333f4lNb4pnL8HxcXFamho0Ny5c4P2Cbfr/02XrmMo1/hKXkvBdItPU3dH2dnZqqioaPf5hSSlpqYqNTXVvz1u3DjdfPPNev311/XMM8+YLvMyGRkZ/n+PGjVKKSkpGjBggDZu3Kj58+d3ej0d8cYbbygjI0Nutzton3C7/uGm29zB9OnTR5GRkaqvrw/YX19fr4SEhDbHJCQkhNS/syxevFjvvfeePvroo5CWoZCknj176tZbb1VlZaWh6kITGxurIUOGBK0nXL8HJ0+e1LZt27RgwYKQxoXb9b90HUO5xlfyWgqm2wRMr169NGbMGJWVlfn3+Xw+lZWVBfwf5r+lpqYG9Jek0tLSoP1NsyxLixcv1jvvvKMPP/xQgwYNCnmO1tZWHTlyRImJiQYqDF1TU5OqqqqC1hNu34NL1qxZo759+2rq1KkhjQu36z9o0CAlJCQEXGOv16t9+/YFvcZX8loKKqRHwmFu/fr1ltPptIqKiqy///3v1kMPPWTFxsZaHo/HsizL+tWvfmU9/vjj/v4ff/yx1aNHD+vFF1+0PvvsMysvL8/q2bOndeTIkS6pf9GiRZbL5bK2b99u1dXV+duXX37p7/PNc1i+fLn1wQcfWFVVVdbBgwetBx54wIqKirKOHj3aFadg/fa3v7W2b99uVVdXWx9//LGVlpZm9enTxzpz5kyb9Yfb98Cy/vMTk/79+1uPPfbYZcfC8fqfP3/eOnz4sHX48GFLkvXyyy9bhw8ftk6ePGlZlmU999xzVmxsrLVlyxbr008/taZPn24NGjTI+te//uWf48c//rH12muv+be/7bX0XXWrgLEsy3rttdes/v37W7169bLGjh1r7d27139swoQJ1pw5cwL6b9y40RoyZIjVq1cva/jw4db777/fyRV/TVKbbc2aNf4+3zyHJUuW+M83Pj7e+ulPf2odOnSo84v/PzNnzrQSExOtXr16Wddff701c+ZMq7Ky0n883L8HlmVZH3zwgSXJOnbs2GXHwvH6f/TRR23+d3OpTp/PZz311FNWfHy85XQ6rUmTJl12bgMGDLDy8vIC9rX3WvquWK4BgDHd5hkMgPBDwAAwhoABYAwBA8AYAgaAMQQMAGMIGADGEDAAjCFgABhDwAAwhoABYAwBA8CY/w8s4xrE6kx3NQAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# The sparse matrix reprenstation of the system\n",
        "plt.imshow(graph.jacobian(ordering)[0][2:,:], cmap=\"RdYlGn\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 212
        },
        "id": "qcekbRZweMlN",
        "outputId": "9a71137a-039e-43a1-f2a9-fd032de0175a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "New Ordering: Position 0: x3, x2, x1, x0, u2, u1, u0\n",
            "\n"
          ]
        },
        {
          "data": {
            "image/svg+xml": [
              "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
              "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
              " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
              "<!-- Generated by graphviz version 2.43.0 (0)\n",
              " -->\n",
              "<!-- Title: %3 Pages: 1 -->\n",
              "<svg width=\"278pt\" height=\"116pt\"\n",
              " viewBox=\"0.00 0.00 278.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
              "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
              "<title>%3</title>\n",
              "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 274,-112 274,4 -4,4\"/>\n",
              "<!-- var8430738502437568512 -->\n",
              "<g id=\"node1\" class=\"node\">\n",
              "<title>var8430738502437568512</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u0</text>\n",
              "</g>\n",
              "<!-- var8430738502437568513 -->\n",
              "<g id=\"node2\" class=\"node\">\n",
              "<title>var8430738502437568513</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u1</text>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&gt;var8430738502437568513 -->\n",
              "<g id=\"edge1\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&gt;var8430738502437568513</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M54.22,-90C56.64,-90 59.11,-90 61.6,-90\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"61.74,-93.5 71.74,-90 61.74,-86.5 61.74,-93.5\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514 -->\n",
              "<g id=\"node3\" class=\"node\">\n",
              "<title>var8430738502437568514</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"171\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">u2</text>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&gt;var8430738502437568514 -->\n",
              "<g id=\"edge2\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&gt;var8430738502437568514</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M54.09,-90C76.47,-90 108.5,-90 133.37,-90\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"133.55,-93.5 143.55,-90 133.55,-86.5 133.55,-93.5\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320 -->\n",
              "<g id=\"node4\" class=\"node\">\n",
              "<title>var8646911284551352320</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&gt;var8646911284551352320 -->\n",
              "<g id=\"edge4\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&gt;var8646911284551352320</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M27,-71.83C27,-64.13 27,-54.97 27,-46.42\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"30.5,-46.41 27,-36.41 23.5,-46.41 30.5,-46.41\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321 -->\n",
              "<g id=\"node5\" class=\"node\">\n",
              "<title>var8646911284551352321</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
              "</g>\n",
              "<!-- var8430738502437568512&#45;&gt;var8646911284551352321 -->\n",
              "<g id=\"edge7\" class=\"edge\">\n",
              "<title>var8430738502437568512&#45;&gt;var8646911284551352321</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M42.27,-74.73C52.2,-64.8 65.32,-51.68 76.44,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"79.16,-42.79 83.75,-33.25 74.21,-37.84 79.16,-42.79\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&gt;var8430738502437568514 -->\n",
              "<g id=\"edge3\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&gt;var8430738502437568514</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M126.22,-90C128.64,-90 131.11,-90 133.6,-90\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"133.74,-93.5 143.74,-90 133.74,-86.5 133.74,-93.5\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&gt;var8646911284551352320 -->\n",
              "<g id=\"edge5\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&gt;var8646911284551352320</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M83.73,-74.73C73.8,-64.8 60.68,-51.68 49.56,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"51.79,-37.84 42.25,-33.25 46.84,-42.79 51.79,-37.84\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&gt;var8646911284551352321 -->\n",
              "<g id=\"edge8\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&gt;var8646911284551352321</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M99,-71.83C99,-64.13 99,-54.97 99,-46.42\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-46.41 99,-36.41 95.5,-46.41 102.5,-46.41\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322 -->\n",
              "<g id=\"node6\" class=\"node\">\n",
              "<title>var8646911284551352322</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"171\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
              "</g>\n",
              "<!-- var8430738502437568513&#45;&gt;var8646911284551352322 -->\n",
              "<g id=\"edge11\" class=\"edge\">\n",
              "<title>var8430738502437568513&#45;&gt;var8646911284551352322</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M114.27,-74.73C124.2,-64.8 137.32,-51.68 148.44,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"151.16,-42.79 155.75,-33.25 146.21,-37.84 151.16,-42.79\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&gt;var8646911284551352320 -->\n",
              "<g id=\"edge6\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&gt;var8646911284551352320</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M149.13,-79.06C124.78,-66.89 85.24,-47.12 57.77,-33.39\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"59.11,-30.14 48.6,-28.8 55.98,-36.4 59.11,-30.14\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&gt;var8646911284551352321 -->\n",
              "<g id=\"edge9\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&gt;var8646911284551352321</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M155.73,-74.73C145.8,-64.8 132.68,-51.68 121.56,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"123.79,-37.84 114.25,-33.25 118.84,-42.79 123.79,-37.84\"/>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&gt;var8646911284551352322 -->\n",
              "<g id=\"edge12\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&gt;var8646911284551352322</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M171,-71.83C171,-64.13 171,-54.97 171,-46.42\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"174.5,-46.41 171,-36.41 167.5,-46.41 174.5,-46.41\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352323 -->\n",
              "<g id=\"node7\" class=\"node\">\n",
              "<title>var8646911284551352323</title>\n",
              "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
              "<text text-anchor=\"middle\" x=\"243\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
              "</g>\n",
              "<!-- var8430738502437568514&#45;&gt;var8646911284551352323 -->\n",
              "<g id=\"edge14\" class=\"edge\">\n",
              "<title>var8430738502437568514&#45;&gt;var8646911284551352323</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M186.27,-74.73C196.2,-64.8 209.32,-51.68 220.44,-40.56\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"223.16,-42.79 227.75,-33.25 218.21,-37.84 223.16,-42.79\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352320&#45;&gt;var8646911284551352321 -->\n",
              "<g id=\"edge10\" class=\"edge\">\n",
              "<title>var8646911284551352320&#45;&gt;var8646911284551352321</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M54.22,-18C56.64,-18 59.11,-18 61.6,-18\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"61.74,-21.5 71.74,-18 61.74,-14.5 61.74,-21.5\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352321&#45;&gt;var8646911284551352322 -->\n",
              "<g id=\"edge13\" class=\"edge\">\n",
              "<title>var8646911284551352321&#45;&gt;var8646911284551352322</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M126.22,-18C128.64,-18 131.11,-18 133.6,-18\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"133.74,-21.5 143.74,-18 133.74,-14.5 133.74,-21.5\"/>\n",
              "</g>\n",
              "<!-- var8646911284551352322&#45;&gt;var8646911284551352323 -->\n",
              "<g id=\"edge15\" class=\"edge\">\n",
              "<title>var8646911284551352322&#45;&gt;var8646911284551352323</title>\n",
              "<path fill=\"none\" stroke=\"black\" d=\"M198.22,-18C200.64,-18 203.11,-18 205.6,-18\"/>\n",
              "<polygon fill=\"black\" stroke=\"black\" points=\"205.74,-21.5 215.74,-18 205.74,-14.5 205.74,-21.5\"/>\n",
              "</g>\n",
              "</g>\n",
              "</svg>\n"
            ],
            "text/plain": [
              "<gtbook.display.show at 0x7e231bfa7a50>"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# The resulting Bayesian graph under a different ordering\n",
        "# The resulting dependency would change, like we have described before\n",
        "\n",
        "# New ordering\n",
        "N = 4\n",
        "ordering = gtsam.Ordering()\n",
        "for k in range(N-1, -1, -1):\n",
        "  ordering.push_back(X(k))\n",
        "for k in range(N-1, -1, -1):\n",
        "  if k > 0:\n",
        "    ordering.push_back(U(k-1))\n",
        "print(\"New Ordering:\", ordering)\n",
        "dag = graph.eliminateSequential(ordering)\n",
        "show(dag, hints={\"u\":1, \"x\":0})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JPGn6rkmolwi"
      },
      "source": [
        "###  Solving the LQR Problem via the Discrete-Time Riccati Equation\n",
        "\n",
        "This code solves a finite-horizon Linear Quadratic Regulator (LQR) problem using the **classical backward Riccati recursion** followed by a **forward rollout** of optimal control inputs.\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "####  Backward Riccati Recursion\n",
        "\n",
        "We compute the optimal feedback gain $K_t$ by solving the **discrete Riccati equation** backward from the terminal cost:\n",
        "\n",
        "$$\n",
        "S_t = R + B^\\top P_{t+1} B\n",
        "$$\n",
        "\n",
        "$$\n",
        "K_t = S_t^{-1} B^\\top P_{t+1} A\n",
        "$$\n",
        "\n",
        "$$\n",
        "P_t = Q + A^\\top P_{t+1} A - A^\\top P_{t+1} B K_t\n",
        "$$\n",
        "\n",
        "The recursion builds the list of gains $K_0, K_1, K_2$, used to compute the optimal control law $u_t = -K_t x_t$.\n",
        "\n",
        "---\n",
        "\n",
        "####  Forward Rollout of the Trajectory\n",
        "\n",
        "Given the initial state $x_0$ and feedback gains:\n",
        "- Compute each optimal control input:  \n",
        "  $u_t = -K_t x_t$\n",
        "- Apply dynamics to get the next state:  \n",
        "  $x_{t+1} = A x_t + B u_t$\n",
        "\n",
        "This yields the full trajectory of states and controls.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 570
        },
        "id": "bZt3HQ5xCVoN",
        "outputId": "b29a789a-92ea-4fff-8087-4e056f85d1b3"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAGGCAYAAACNL1mYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbRRJREFUeJzt3XlcVPX+x/HXsMyArC6AG4lLmbuGu5lalHtZlra6V5paZnWvdkuzzWtdy36pWZba5tXUMlPLlNSy676UlpmZa4poKggoy8z5/TEyMMygjArD8n4+HueBnPmecz4DJ+LN93u+X5NhGAYiIiIiIiIiUih8vF2AiIiIiIiISGmm4C0iIiIiIiJSiBS8RURERERERAqRgreIiIiIiIhIIVLwFhERERERESlECt4iIiIiIiIihUjBW0RERERERKQQKXiLiIiIiIiIFCIFbxEREREREZFCpOAtIuJlHTt2pGPHjt4uo0x54YUXMJlM3i7jinnzfcTExDBgwACvXLu0OXDgACaTiTlz5ni7lMsyYMAAYmJivF2GiEixpuAtInKVzZkzB5PJ5Nj8/PyoVq0aAwYM4K+//vJ2eW5Nnz69xP7Sn5+0tDReeOEF1qxZ4+1SyrSUlBTGjx9Pw4YNCQoKomLFijRt2pQnnniCo0ePOtotX76cF1544Yqu9eqrr7J48eIrKziX7D9sXGor7n84O3r0KC+88AI7duzwdikiImWWn7cLEBEprV588UVq1qzJ+fPn2bBhA3PmzGHdunXs2rWLgIAAR7tvv/3Wi1XaTZ8+nUqVKpWqHsy0tDQmTJgA4BKMnnvuOcaMGeOFqsqWzMxMbrrpJn777Tf69+/PyJEjSUlJ4ZdffmHu3LnceeedVK1aFbAH72nTpl1R+H711Ve5++676dWr11Wp/6677qJOnTqOz1NSUhg2bBh33nknd911l2N/VFTUFV2nRo0anDt3Dn9//ys6T36OHj3KhAkTiImJoWnTplf9/DNnzsRms13184qIlCYK3iIihaRr1640b94cgCFDhlCpUiUmTZrEkiVL6NOnj6Od2Wz2Volllp+fH35++l9gYVu8eDHbt2/n008/5f7773d67fz582RkZHipsoJp3LgxjRs3dnx+8uRJhg0bRuPGjXnwwQfzPe78+fOYzWZ8fAo2sNBkMjn9Ma6kSE1NJSgoqND+YCAiUppoqLmISBFp3749APv27XPa7+4Z7/Pnz/PCCy9w3XXXERAQQJUqVbjrrrucjrXZbLz11ls0atSIgIAAIiIi6NKlC1u2bHG0mT17NjfffDORkZFYLBbq16/PO++843StmJgYfvnlF9auXevR0Nm//vqLQYMGERUVhcVioUGDBsyaNcul3ZEjR+jVqxdBQUFERkby5JNPsmLFCkwmk9Mw8PyeGc779cnIyGDcuHHExsYSFhZGUFAQ7du3Z/Xq1Y42Bw4cICIiAoAJEyY43ld2b6q7Z6OzsrJ46aWXqF27NhaLhZiYGJ599lnS09Ndvl49evRg3bp1tGzZkoCAAGrVqsVHH310ya8ZwNmzZxk1ahQxMTFYLBYiIyO59dZb2bZtm1O7BQsWEBsbS2BgIJUqVeLBBx+85KMKDRs2pFOnTi77bTYb1apV4+6773baN2XKFBo0aEBAQABRUVE8+uijnD592ulYwzB4+eWXqV69OuXKlaNTp0788ssvBXqv2fdru3btXF4LCAggNDQUsD8jPG3aNACnIdzZ/vOf/9C2bVsqVqxIYGAgsbGxLFy40Ol8JpOJ1NRUPvzwQ8fxue+ngt6vnlqzZg0mk4l58+bx3HPPUa1aNcqVK0dycjKnTp3i6aefplGjRgQHBxMaGkrXrl356aefnM6R3zPev/32G3fffTcVKlQgICCA5s2bs2TJEpcazpw5w5NPPum4p6pXr06/fv04efIka9asoUWLFgAMHDjQ8bXJfa2C3GsDBgwgODiYffv20a1bN0JCQnjggQccr+V9xrug99eWLVvo3LkzlSpVIjAwkJo1azJo0CBPvgUiIiWC/twvIlJEDhw4AED58uUv2s5qtdKjRw/i4+O59957eeKJJzh79iwrV65k165d1K5dG4DBgwczZ84cunbtypAhQ8jKyuKHH35gw4YNjp72d955hwYNGnD77bfj5+fHV199xWOPPYbNZmP48OEATJkyhZEjRxIcHMy//vUv4NJDZ48fP07r1q0xmUyMGDGCiIgIvv76awYPHkxycjKjRo0C4Ny5c9xyyy0cOnSIxx9/nKpVq/Lxxx/z3XffXe6XkeTkZN5//33uu+8+Hn74Yc6ePcsHH3xA586d2bRpE02bNiUiIoJ33nnHZVhw7t7LvIYMGcKHH37I3XffzVNPPcXGjRuZOHEiu3fv5osvvnBq+8cff3D33XczePBg+vfvz6xZsxgwYACxsbE0aNDgovUPHTqUhQsXMmLECOrXr8/ff//NunXr2L17NzfccANgnydg4MCBtGjRgokTJ3L8+HHeeustfvzxR7Zv3054eLjbc/ft25cXXniBhIQEKleu7Ni/bt06jh49yr333uvY9+ijjzqu8/jjj7N//36mTp3K9u3b+fHHHx29mOPGjePll1+mW7dudOvWjW3btnHbbbcVqLe6Ro0aAHz00Uc899xz+U4E9+ijj3L06FFWrlzJxx9/7PL6W2+9xe23384DDzxARkYG8+bN45577mHp0qV0794dgI8//pghQ4bQsmVLHnnkEQDHfysFvV+vxEsvvYTZbObpp58mPT0ds9nMr7/+yuLFi7nnnnuoWbMmx48f591336VDhw78+uuvjmH27vzyyy+0a9eOatWqMWbMGIKCgvjss8/o1asXixYt4s477wTsw9/bt2/P7t27GTRoEDfccAMnT55kyZIlHDlyhHr16vHiiy8ybtw4HnnkEccfANu2bQt4dq9lZWXRuXNnbrzxRv7zn/9Qrly5fOsvyP2VmJjIbbfdRkREBGPGjCE8PJwDBw7w+eefX/H3Q0Sk2DFEROSqmj17tgEYq1atMk6cOGEcPnzYWLhwoREREWFYLBbj8OHDTu07dOhgdOjQwfH5rFmzDMB44403XM5ts9kMwzCM7777zgCMxx9/PN82hmEYaWlpLq937tzZqFWrltO+Bg0aONVwKYMHDzaqVKlinDx50mn/vffea4SFhTmuO2XKFAMwPvvsM0eb1NRUo06dOgZgrF692rG/Ro0aRv/+/V2ulffrk5WVZaSnpzu1OX36tBEVFWUMGjTIse/EiRMGYIwfP97lnOPHjzdy/y9wx44dBmAMGTLEqd3TTz9tAMZ3333nVCdgfP/99459iYmJhsViMZ566imXa+UVFhZmDB8+PN/XMzIyjMjISKNhw4bGuXPnHPuXLl1qAMa4cePyfR979uwxAOPtt992Oudjjz1mBAcHO74vP/zwgwEYn376qVO7b775xml/YmKiYTabje7duzvdV88++6wBuP1+5ZaWlmbUrVvXAIwaNWoYAwYMMD744APj+PHjLm2HDx9u5PdrSd77OCMjw2jYsKFx8803O+0PCgpyW1NB79dLcXdPrV692gCMWrVquZzn/PnzhtVqddq3f/9+w2KxGC+++KLTPsCYPXu2Y98tt9xiNGrUyDh//rxjn81mM9q2bWtce+21jn3jxo0zAOPzzz93qTf7e7Z582aX8xuGZ/da//79DcAYM2aMy3X69+9v1KhRw/F5Qe+vL774wgCMzZs3u5xTRKS00VBzEZFCEhcXR0REBNHR0dx9990EBQWxZMkSqlevftHjFi1aRKVKlRg5cqTLa9k9hosWLcJkMjF+/Ph82wAEBgY6/p2UlMTJkyfp0KEDf/75J0lJSZf1vgzDYNGiRfTs2RPDMDh58qRj69y5M0lJSY5h08uXL6dKlSpOQ5zLlSvn6JG8HL6+vo7n4m02G6dOnSIrK4vmzZu7DNcuqOXLlwMwevRop/1PPfUUAMuWLXPaX79+fUfPIUBERAR169blzz//vOS1wsPD2bhxo9OM3rlt2bKFxMREHnvsMafnfrt3787111/vUktu1113HU2bNmX+/PmOfVarlYULF9KzZ0/H/bBgwQLCwsK49dZbnb5/sbGxBAcHO4btr1q1ioyMDEaOHOl0XxW0hzgwMJCNGzfyzDPPAPbe1cGDB1OlShVGjhzpMoz/YufJdvr0aZKSkmjfvn2Bvt+e3K9Xon///k51AlgsFsdz3larlb///pvg4GDq1q170WueOnWK7777jj59+nD27FlHvX///TedO3dm7969jqHgixYtokmTJo4e8NwutdTc5dxrw4YNu+g5oeD3V3Zv+tKlS8nMzLzkeUVESjINNRcRKSTTpk3juuuuIykpiVmzZvH9999jsVguedy+ffuoW7fuRSf/2rdvH1WrVqVChQoXPdePP/7I+PHjWb9+PWlpaU6vJSUlERYWlu+xVquVEydOOO2rUKECZ86c4cyZM7z33nu89957bo9NTEwE4ODBg9SpU8clANStW/eidV/Khx9+yOTJk/ntt9+cfmGvWbPmZZ3v4MGD+Pj4OM1gDVC5cmXCw8M5ePCg0/5rrrnG5Rzly5d3PL+a39fObDbz2muv0b9/f6Kjo4mNjaVbt27069ePWrVqOWoB91+j66+/nnXr1l30vfTt25dnn32Wv/76i2rVqrFmzRoSExPp27evo83evXtJSkoiMjLS7Tlyf/8Arr32WqfXIyIiLvnIRLawsDBee+01XnvtNQ4ePEh8fDz/+c9/mDp1KmFhYbz88suXPMfSpUt5+eWX2bFjh1NYL8ga5idOnCjw/Xol3N172fMwTJ8+nf3792O1Wh2vVaxYMd9z/fHHHxiGwfPPP8/zzz+fb83VqlVj37599O7d+7Jq9vRe8/Pzu+QfDqHg91eHDh3o3bs3EyZM4M0336Rjx4706tWL+++/v0A/K0VEShIFbxGRQtKyZUvHs9a9evXixhtv5P7772fPnj0EBwcX+vX37dvHLbfcwvXXX88bb7xBdHQ0ZrOZ5cuX8+abb15y+Z/Dhw+7hInVq1dz/fXXA/Dggw/Sv39/t8de7Fnq/OQXoqxWK76+vo7PP/nkEwYMGECvXr145plniIyMxNfXl4kTJ7pMXHe1asgrdz25GYYB5P+169ixI3369KF9+/Z88cUXfPvtt7z++utMmjSJzz//nK5du15R/WAP3mPHjmXBggWMGjWKzz77jLCwMLp06eJoY7PZiIyM5NNPP3V7juyJ6a62GjVqMGjQIO68805q1arFp59+esng/cMPP3D77bdz0003MX36dKpUqYK/vz+zZ89m7ty5l7xm9n1+te/XvPL2doN9ebPnn3+eQYMG8dJLL1GhQgV8fHwYNWrURf/7y37t6aefpnPnzm7b5P0jUVHI3YN/MQW9v0wmEwsXLmTDhg189dVXrFixgkGDBjF58mQ2bNhQJD8nRUSKioK3iEgRyA6GnTp1YurUqRddQ7p27dps3LiRzMzMfJfpqV27NitWrODUqVP59np/9dVXpKens2TJEqce2tyzf2dzFzgrV67MypUrnfY1adKE0NBQQkJCsFqtxMXF5fs+wB60du3ahWEYTtfYs2ePS9vy5ctz5swZl/0HDx509AYDLFy4kFq1avH55587nTPvsPuChujsOm02G3v37qVevXqO/cePH+fMmTOOScIKKr+vXbYqVarw2GOP8dhjj5GYmMgNN9zAK6+8QteuXR3X2rNnDzfffLPTOfbs2XPJWmrWrEnLli2ZP38+I0aM4PPPP6dXr15OPYi1a9dm1apVtGvXzm1gzJZ9rb179zp9D06cOOEyO7UnypcvT+3atdm1a5djX37fr0WLFhEQEMCKFSuc3sPs2bNd2ro7R0RERIHv16tt4cKFdOrUiQ8++MBp/5kzZ6hUqVK+x2V/rf39/S9Zc96vozv5fW2v9F67WE0Fub+ytW7dmtatW/PKK68wd+5cHnjgAebNm8eQIUMu6/oiIsWRnvEWESkiHTt2pGXLlkyZMoXz58/n2653796cPHmSqVOnuryW3aPau3dvDMNgwoQJ+bbJ7pXN/hzsw8vdBZagoCCX0BsQEEBcXJzTVr58eXx9fenduzeLFi1y+wt/7iHW3bp14+jRo05LP6Wlpbkd8lu7dm02bNjgNFv20qVLOXz4sFM7d+9r48aNrF+/3qld9ozL7sJ8Xt26dQPsM7zn9sYbbwA4Zs4uqPy+dlar1eXZ+sjISKpWreoYQt28eXMiIyOZMWOG07Dqr7/+mt27dxeolr59+7JhwwZmzZrFyZMnnYaZA/Tp0wer1cpLL73kcmxWVpbjaxYXF4e/vz9vv/2209c779cpPz/99BMnT5502X/w4EF+/fVXpyHOQUFBgOv3y9fXF5PJ5DRM+8CBAyxevNjlvO7uY0/u16vN19fX6esG9uefL7UsXGRkJB07duTdd9/l2LFjLq/nrrl379789NNPLjPvQ85/I/l9ba/GveZOQe+v06dPu3x9mjZtClDg5/9FREoK9XiLiBShZ555hnvuuYc5c+YwdOhQt2369evHRx99xOjRo9m0aRPt27cnNTWVVatW8dhjj3HHHXfQqVMnHnroIf7v//6PvXv30qVLF2w2Gz/88AOdOnVixIgR3HbbbZjNZnr27Mmjjz5KSkoKM2fOJDIy0uWX+djYWN555x1efvll6tSpQ2RkpEsPWG7//ve/Wb16Na1ateLhhx+mfv36nDp1im3btrFq1SpOnToFwMMPP8zUqVPp168fW7dupUqVKnz88cdulyEaMmQICxcupEuXLvTp04d9+/bxySefOJaEytajRw8+//xz7rzzTrp3787+/fuZMWMG9evXJyUlxdEuMDCQ+vXrM3/+fK677joqVKhAw4YNadiwocu1mzRpQv/+/Xnvvfc4c+YMHTp0YNOmTXz44Yf06tXL7drYl+Ps2bNUr16du+++myZNmhAcHMyqVavYvHkzkydPBuy9nJMmTWLgwIF06NCB++67z7HEU0xMDE8++eQlr9OnTx+efvppnn76aSpUqODSa9qhQwceffRRJk6cyI4dO7jtttvw9/dn7969LFiwgLfeeou7776biIgInn76aSZOnEiPHj3o1q0b27dv5+uvv75oj222lStXMn78eG6//XZat25NcHAwf/75J7NmzSI9Pd2xrjrY70GAxx9/nM6dO+Pr68u9995L9+7deeONN+jSpQv3338/iYmJTJs2jTp16vDzzz87XS82NpZVq1bxxhtvULVqVWrWrEmrVq0KfL9ebT169ODFF19k4MCBtG3blp07d/Lpp586jR7Iz7Rp07jxxhtp1KgRDz/8MLVq1eL48eOsX7+eI0eOONYCf+aZZ1i4cCH33HMPgwYNIjY2llOnTrFkyRJmzJhBkyZNqF27NuHh4cyYMYOQkBCCgoJo1aoVNWvWvOJ7zZ2C3l8ffvgh06dP584776R27dqcPXuWmTNnEhoa6vhjmIhIqVH0E6mLiJRu2cuJuVsix2q1GrVr1zZq165tZGVlGYbhulyWYdiXT/rXv/5l1KxZ0/D39zcqV65s3H333ca+ffscbbKysozXX3/duP766w2z2WxEREQYXbt2NbZu3epos2TJEqNx48ZGQECAERMTY0yaNMmxXNn+/fsd7RISEozu3bsbISEhBlCgpcWOHz9uDB8+3IiOjnbUeMsttxjvvfeeU7uDBw8at99+u1GuXDmjUqVKxhNPPOFYVij3cmKGYRiTJ082qlWrZlgsFqNdu3bGli1bXL4+NpvNePXVV40aNWoYFovFaNasmbF06VKXJY0MwzD+97//GbGxsYbZbHZaBirvMlyGYRiZmZnGhAkTHF/z6OhoY+zYsU7LORmGfTmx7t27u3w93H0f80pPTzeeeeYZo0mTJkZISIgRFBRkNGnSxJg+fbpL2/nz5xvNmjUzLBaLUaFCBeOBBx4wjhw54tTG3fvI1q5dO7dLpOX23nvvGbGxsUZgYKAREhJiNGrUyPjHP/5hHD161NHGarUaEyZMMKpUqWIEBgYaHTt2NHbt2pXv8m+5/fnnn8a4ceOM1q1bG5GRkYafn58RERFhdO/e3WmJNsOw388jR440IiIiDJPJ5PS+PvjgA+Paa681LBaLcf311xuzZ892+95/++0346abbjICAwNdljsr6P16MRdbTmzBggUu7c+fP2889dRTjq9du3btjPXr17vcK+6WEzMMw9i3b5/Rr18/o3Llyoa/v79RrVo1o0ePHsbChQud2v3999/GiBEjjGrVqhlms9moXr260b9/f6fl07788kujfv36hp+fn8u1CnKv9e/f3wgKCnL7dXH3355hXPr+2rZtm3HfffcZ11xzjWGxWIzIyEijR48expYtW9xeR0SkJDMZRp4xPiIiIoVszZo1dOrUyTHhmEhZtm/fPurUqcPHH3/Mgw8+6O1yRESkEOgZbxEREREvyn70oyDD90VEpGTSM94iIiIiXjJr1ixmzZpFuXLlaN26tbfLERGRQqIebxEREREveeSRRzh16hQLFiwgPDzc2+WIiEgh0TPeIiIiIiIiIoVIPd4iIiIiIiIihUjBW0RERERERKQQlYjJ1Ww2G0ePHiUkJASTyeTtckRERERERKSMMwyDs2fPUrVqVXx8Lt6nXSKC99GjR4mOjvZ2GSIiIiIiIiJODh8+TPXq1S/apkQE75CQEMD+hkJDQ71cjYiIiIiIiJR1ycnJREdHO/LqxZSI4J09vDw0NFTBW0RERERERIqNgjwOrcnVRERERERERAqRgreIiIiIiIhIIVLwFhERERERESlEJeIZbxERERERkdLMZrORkZHh7TIkF39/f3x9fa/KuRS8RUREREREvCgjI4P9+/djs9m8XYrkER4eTuXKlQs0gdrFKHiLiIiIiIh4iWEYHDt2DF9fX6Kjo/Hx0dPAxYFhGKSlpZGYmAhAlSpVruh8Ct4iIiIiIiJekpWVRVpaGlWrVqVcuXLeLkdyCQwMBCAxMZHIyMgrGnbu8Z9Tvv/+e3r27EnVqlUxmUwsXrz4ksesWbOGG264AYvFQp06dZgzZ85llCoiIiIiIlK6WK1WAMxms5crEXey/xiSmZl5RefxOHinpqbSpEkTpk2bVqD2+/fvp3v37nTq1IkdO3YwatQohgwZwooVKzwuVkREREREpDS60meIpXBcre+Lx0PNu3btSteuXQvcfsaMGdSsWZPJkycDUK9ePdatW8ebb75J586dPb188WbNgowUCAz3diUiIiIiIiJSTBT6k/vr168nLi7OaV/nzp1Zv359vsekp6eTnJzstJUIG2fA1Obw0zwwDG9XIyIiIiIiUqKsWbMGk8nEmTNnLtouJiaGKVOmFElNV0OhB++EhASioqKc9kVFRZGcnMy5c+fcHjNx4kTCwsIcW3R0dGGXeeVsNti1CFJPwBePwpwekPibt6sSERERERG56gYMGIDJZMJkMmE2m6lTpw4vvvgiWVlZV3Tetm3bcuzYMcLCwgCYM2cO4eHhLu02b97MI488ckXXKkrFcq76sWPHkpSU5NgOHz7s7ZIuzccHBq2AW8aBXyAcXAcz2sGqFyAjzdvViYiIiIiIXFVdunTh2LFj7N27l6eeeooXXniB119//YrOaTabC7RudkRERImaBb7Qg3flypU5fvy4077jx48TGhrqmJ49L4vFQmhoqNNWIviZof1TMHwjXNcVbFmw7k2Y1gr2fO3t6kRERERERK4ai8VC5cqVqVGjBsOGDSMuLo4lS5Zw+vRp+vXrR/ny5SlXrhxdu3Zl7969juMOHjxIz549KV++PEFBQTRo0IDly5cDzkPN16xZw8CBA0lKSnL0rr/wwguA61DzQ4cOcccddxAcHExoaCh9+vRxyqEvvPACTZs25eOPPyYmJoawsDDuvfdezp49WyRfq0IP3m3atCE+Pt5p38qVK2nTpk1hX9p7yteA++fBvXMhLBqSDsF/74X/3gdnDnm7OhERERERKaYMwyAtI8srm3GF81QFBgaSkZHBgAED2LJlC0uWLGH9+vUYhkG3bt0cS3INHz6c9PR0vv/+e3bu3MmkSZMIDg52OV/btm2ZMmUKoaGhHDt2jGPHjvH000+7tLPZbNxxxx2cOnWKtWvXsnLlSv7880/69u3r1G7fvn0sXryYpUuXsnTpUtauXcu///3vK3rPBeXxrOYpKSn88ccfjs/379/Pjh07qFChAtdccw1jx47lr7/+4qOPPgJg6NChTJ06lX/84x8MGjSI7777js8++4xly5ZdvXdRXF3fHWp1hLWvwfqpsGc57FsNHf4BbUbYe8hFREREREQuOJdppf447yy9/OuLnSln9jgiYhgG8fHxrFixgq5du7J48WJ+/PFH2rZtC8Cnn35KdHQ0ixcv5p577uHQoUP07t2bRo0aAVCrVi235zWbzYSFhWEymahcuXK+14+Pj2fnzp3s37/fMT/YRx99RIMGDdi8eTMtWrQA7AF9zpw5hISEAPDQQw8RHx/PK6+84vF79pTHPd5btmyhWbNmNGvWDIDRo0fTrFkzxo0bB8CxY8c4dCinV7dmzZosW7aMlStX0qRJEyZPnsz7779f+pYSy485CG6dAEPXQY12kHUO4ifAu+3hwDpvVyciIiIiInJZli5dSnBwMAEBAXTt2pW+ffsyYMAA/Pz8aNWqlaNdxYoVqVu3Lrt37wbg8ccf5+WXX6Zdu3aMHz+en3/++Yrq2L17N9HR0U6TctevX5/w8HDHNcE+PD07dANUqVKFxMTEK7p2QXn854yOHTtedAjCnDlz3B6zfft2Ty9VukTWgwHL7EuNffscnPgN5nSHxvfCbS9DcIS3KxQRERERES8L9Pfl1xe900kZ6O/rUftOnTrxzjvvYDabqVq1Kn5+fixZsuSSxw0ZMoTOnTuzbNkyvv32WyZOnMjkyZMZOXLk5ZZeIP7+/k6fm0wmbDZboV4zW7Gc1bzUMpmg6X0wYjM0HwSY4Od5MDUWNn8ANqu3KxQRERERES8ymUyUM/t5ZbvUTOJ5BQUFUadOHa655hr8/Ox9uvXq1SMrK4uNGzc62v3999/s2bOH+vXrO/ZFR0czdOhQPv/8c5566ilmzpzp9hpmsxmr9eI5qV69ehw+fNhpNaxff/2VM2fOOF3TmxS8vaFcBejxJgxZBZUbw/kkWDYaPrgVju7wdnUiIiIiIiKX5dprr+WOO+7g4YcfZt26dfz00088+OCDVKtWjTvuuAOAUaNGsWLFCvbv38+2bdtYvXo19erVc3u+mJgYUlJSiI+P5+TJk6SluS7VHBcXR6NGjXjggQfYtm0bmzZtol+/fnTo0IHmzZsX6vstKAVvb6reHB5eDV0mgTkE/toKMzvB8mfsYVxERERERKSEmT17NrGxsfTo0YM2bdpgGAbLly93DPW2Wq0MHz6cevXq0aVLF6677jqmT5/u9lxt27Zl6NCh9O3bl4iICF577TWXNiaTiS+//JLy5ctz0003ERcXR61atZg/f36hvk9PmIwrnTO+CCQnJxMWFkZSUlLJWdPbU2cTYMWzsGuR/fPgKOj8KjTsbR+iLiIiIiIipc758+fZv38/NWvWJCAgwNvlSB4X+/54klPV411chFSGu2fBQ4uhYh1IOQ6LBsNHd8DJvZc8XERERERERIonBe/ipnYnGPY/6PQc+AXA/rXwTlv47mXIPOft6kRERERERMRDCt7FkZ8FOjwDj22AOreCNQO+fx2mtYLfv/V2dSIiIiIiIuIBBe/irEJNeGAB9PkIQqrCmYMw9x6Y/yAkHfF2dSIiIiIiIlIACt7FnckE9e+AEZugzQgw+cLur2BqS/jx/8Ca6e0KRURERERE5CIUvEsKSwh0fgWG/gDRrSEzFVY+D+/eBAfXe7s6ERERERERyYeCd0kT1QAGfg23T4XACpD4K8zuAouHQ+pJb1cnIiIiIiIieSh4l0Q+PnDDQzByK9zQz75vxycwtTlsnQM2m1fLExERERERkRwK3iVZuQpw+9sw6FuIagjnTsNXT8CszpCw09vViYiIiIiICArepcM1reCRtdD5VTAHw5FN9me/vxkL6We9XZ2IiIiIiIiTmJgYpkyZUmzPd7UpeJcWvn7QZjiM2Az1e4Fhgw3TYWoL+OULMAxvVygiIiIiIqVAz5496dKli9vXfvjhB0wmEz///HOR1rR582YeeeQRx+cmk4nFixcXaQ0Xo+Bd2oRWhT4fwgOLoHxNOHsMFgyAT3rD3/u8XZ2IiIiIiJRwgwcPZuXKlRw5csTltdmzZ9O8eXMaN25cpDVFRERQrly5Ir2mJxS8S6tr4+Cx9dBhDPiaYV88TG8DqydC5nlvVyciIiIiIiVUjx49iIiIYM6cOU77U1JSWLBgAYMHD2bdunW0b9+ewMBAoqOjefzxx0lNTc33nIcOHeKOO+4gODiY0NBQ+vTpw/Hjx53afPXVV7Ro0YKAgAAqVarEnXfe6Xgt91DzmJgYAO68805MJhMxMTEcOHAAHx8ftmzZ4nTOKVOmUKNGDWyFPEG1gndp5h8IncbCYxugViewpsPaf8M7beCPeG9XJyIiIiIieRkGZKR6Zyvg46l+fn7069ePOXPmYOQ6ZsGCBVitVtq0aUOXLl3o3bs3P//8M/Pnz2fdunWMGDHC7flsNht33HEHp06dYu3ataxcuZI///yTvn37OtosW7aMO++8k27durF9+3bi4+Np2bKl2/Nt3rwZsPe+Hzt2jM2bNxMTE0NcXByzZ892ajt79mwGDBiAj0/hRmOTYRT/h3+Tk5MJCwsjKSmJ0NBQb5dTMhmG/Vnvb8ZCSoJ9X4M77ROyhVb1bm0iIiIiImXU+fPn2b9/PzVr1iQgIMAegF/10u/nzx4Fc1CBmv7222/Uq1eP1atX07FjRwBuuukmatSogcViwdfXl3fffdfRft26dXTo0IHU1FQCAgKIiYlh1KhRjBo1ipUrV9K1a1f2799PdHQ0AL/++isNGjRg06ZNtGjRgrZt21KrVi0++eQTt/XkPh/Yn/H+4osv6NWrl6PNZ599xtChQzl27BgWi4Vt27bRvHlz/vzzT0cveV4u359cPMmp6vEuK0wmaHiXffK11o+ByccexKe2gPXTwZrl7QpFRERERKSEuP7662nbti2zZs0C4I8//uCHH35g8ODB/PTTT8yZM4fg4GDH1rlzZ2w2G/v373c51+7du4mOjnaEboD69esTHh7O7t27AdixYwe33HLLFdXcq1cvfH19+eKLLwCYM2cOnTp1yjd0X01+hX4FKV4CQqHLRGhyHywbDUc2w4qxsGMu9HgDot0P1xARERERkSLgX87e8+yta3tg8ODBjBw5kmnTpjF79mxq165Nhw4dSElJ4dFHH+Xxxx93Oeaaa665rNICAwMv67jczGYz/fr1Y/bs2dx1113MnTuXt95664rPWxAK3mVVlcYw6FvY9iGsegGO74QPboUb+kHcBChXwdsVioiIiIiUPSZTgYd7e1ufPn144oknmDt3Lh999BHDhg3DZDJxww038Ouvv1KnTp0CnadevXocPnyYw4cPOw01P3PmDPXr1wegcePGxMfHM3DgwAKd09/fH6vV6rJ/yJAhNGzYkOnTp5OVlcVdd91VwHd7ZTTUvCzz8YHmA2HkVmj6gH3fto9ganPY/gkU8sx+IiIiIiJScgUHB9O3b1/Gjh3LsWPHGDBgAAD//Oc/+d///seIESPYsWMHe/fu5csvv8x3crW4uDgaNWrEAw88wLZt29i0aRP9+vWjQ4cONG/eHIDx48fz3//+l/Hjx7N792527tzJpEmT8q0tJiaG+Ph4EhISOH36tGN/vXr1aN26Nf/85z+57777rkpPekEoeAsEVYJe02Hg1xBRD9L+hi+Hw5xucPwXb1cnIiIiIiLF1ODBgzl9+jSdO3emalX7pHCNGzdm7dq1/P7777Rv355mzZoxbtw4x+t5mUwmvvzyS8qXL89NN91EXFwctWrVYv78+Y42HTt2ZMGCBSxZsoSmTZty8803s2nTpnzrmjx5MitXriQ6OppmzZq51JyRkcGgQYOuwlegYDSruTizZsKGd2DNvyEzFUy+0OYx+3rglmBvVyciIiIiUqpcbNZsKRwvvfQSCxYs4Oeff75kW81qLoXD1x/aPQ4jNsH1PcCwwv/ehmkt4dclBV7bT0REREREpDhJSUlh165dTJ06lZEjRxbptRW8xb2w6nDvp3D/ZxB+DST/BZ89BHP7wCnXJQBERERERESKsxEjRhAbG0vHjh2LdJg5KHjLpVzXGR7bCO2fBh9/2PstTG8Na1+HrHRvVyciIiIiIlIgc+bMIT09nfnz5+Pr61uk11bwlkszl4NbnofH1kPNmyDrPKx+Gd5pC3+u8XZ1IiIiIiIixZqCtxRcpWuh3xK4630IioS//4CP7oCFg+FsgrerExERERERKZYUvMUzJhM0vgdGbIaWj4DJB3YthKktYOO7YHNdpF5ERERERC6uBCw2VSbZbLarch4tJyZX5uh2WPqk/SNAlSbQ/U2oHuvdukRERERESgCr1crevXspV64cERERmEwmb5ck2P8QkpGRwYkTJ7BarVx77bX4+Dj3W3uSUxW85crZrLB1Nqx6EdKTABM0Hwi3jIPA8t6uTkRERESkWEtJSeHIkSPq9S6GypUrR5UqVTCbzS6vKXiLd6QkwrfPw8/z7J8HRcBtL0PjvvYh6iIiIiIi4pbVaiUzM9PbZUguvr6++Pn55TsKQcFbvGv/D7DsKTi5x/55jRuh+2SIvN67dYmIiIiIiFwlnuRUTa4mV1/N9jB0HdwyHvwC4eA6mNEOVo6HjFRvVyciIiIiIlKkFLylcPiZof1oGL4RrusKtiz4cQpMaw2/Lfd2dSIiIiIiIkVGwVsKV/kacP88uPe/EBYNSYdg3n3w3/vg9EFvVyciIiIiIlLoFLylaFzfzd77feOT4OMHe5bDtFbwwxuQleHt6kRERERERAqNgrcUHXMQxL0AQ3+0T7iWdQ7iJ8CMG+0TsomIiIiIiJRCCt5S9CKvhwFLodcMKFfJPvv5hz3g80fsS5KJiIiIiIiUIgre4h0mEzS9D0ZugeaDABP8PB+mNofN74PN6u0KRURERERErgoFb/GuwPLQ400YEg9VmsD5JPsa4O/HwdHt3q5ORERERETkiil4S/FQPRYeXg1dXwdLKBzdBjNvhuXP2MO4iIiIiIhICXVZwXvatGnExMQQEBBAq1at2LRp00XbT5kyhbp16xIYGEh0dDRPPvkk58+fv6yCpRTz8YVWj8CIzdDwbjBssOk9eLs5/LwADMPbFYqIiIiIiHjM4+A9f/58Ro8ezfjx49m2bRtNmjShc+fOJCa6nxRr7ty5jBkzhvHjx7N7924++OAD5s+fz7PPPnvFxUspFVIZ7v4AHloMFetAaiJ8PgQ+ugNO7vV2dSIiIiIiIh4xGYZn3YitWrWiRYsWTJ06FQCbzUZ0dDQjR45kzJgxLu1HjBjB7t27iY+Pd+x76qmn2LhxI+vWrSvQNZOTkwkLCyMpKYnQ0FBPypWSLisdfvw/+OE/kHUefPyh3RNw09PgH+jt6kREREREpIzyJKd61OOdkZHB1q1biYuLyzmBjw9xcXGsX7/e7TFt27Zl69atjuHof/75J8uXL6dbt26eXFrKKj8LdHgGHtsA194Gtkx7CJ/WCn7/1tvViYiIiIiIXJKfJ41PnjyJ1WolKirKaX9UVBS//fab22Puv/9+Tp48yY033ohhGGRlZTF06NCLDjVPT08nPT3d8XlycrInZUppVKEm3P8Z7P4KvhkDZw7C3Hvg+h7QdRKEVfd2hSIiIiIiIm4V+qzma9as4dVXX2X69Ols27aNzz//nGXLlvHSSy/le8zEiRMJCwtzbNHR0YVdppQEJhPUvx2Gb4K2I8HkC78thakt4ce3wJrp7QpFRERERERcePSMd0ZGBuXKlWPhwoX06tXLsb9///6cOXOGL7/80uWY9u3b07p1a15//XXHvk8++YRHHnmElJQUfHxcs7+7Hu/o6Gg94y3Ojv8CS0fD4Q32zyPrQ/c3oEYb79YlIiIiIiKlXqE94202m4mNjXWaKM1msxEfH0+bNu7DTlpamku49vX1BSC/zG+xWAgNDXXaRFxENYCBX8Md0yCwAiT+CrO7wOLHIPWkt6sTEREREREBLmOo+ejRo5k5cyYffvghu3fvZtiwYaSmpjJw4EAA+vXrx9ixYx3te/bsyTvvvMO8efPYv38/K1eu5Pnnn6dnz56OAC5y2Xx8oNmDMHIr3NDfvm/Hp/B2LGydAzabV8sTERERERHxaHI1gL59+3LixAnGjRtHQkICTZs25ZtvvnFMuHbo0CGnHu7nnnsOk8nEc889x19//UVERAQ9e/bklVdeuXrvQqRcBbj9/+whfOloOL4TvnoCtn9iH35epbG3KxQRERERkTLK43W8vUHreItHrFmw6T1Y/QpkpIDJB1o+Cp2ehQDdPyIiIiIicuUK7RlvkRLB1w/aPAYjNkP9XmDYYOM7MK0l7Pociv/fmkREREREpBRR8JbSK7Qq9PkQHlwE5WvC2WOwcCB8chf8vc/b1YmIiIiISBmh4C2lX504eGwDdBwLvhbY9x1MbwOrJ0LmeW9XJyIiIiIipZyCt5QN/gHQcQw8th5q3wzWdFj7b5jeGv5Y5e3qRERERESkFFPwlrKlYm148HO4ezYEV4bT++GT3vBZf0g+6u3qRERERESkFFLwlrLHZIKGd9knX2v9mH3W818Xw9QWsH6afVZ0ERERERGRq0TBW8qugFDoMhEeWQvVW9iXHlvxLLzXAQ5t9HZ1IiIiIiJSSih4i1RpDIO+hZ5vQUA4HN8Fs26DJSMh7ZS3qxMRERERkRJOwVsEwMcHYgfAyK3Q9EH7vm0fwduxsO1jsNm8Wp6IiIiIiJRcCt4iuQVVgl7TYOA3EFkfzp2CJSNgdlc4/ou3qxMRERERkRJIwVvEnRpt4NHv4daXwD8IDm+AGe1hxb8gPcXb1YmIiIiISAmi4C2SH19/aPc4jNgE9XqCYYX1U2FaS/h1CRiGtysUEREREZESQMFb5FLCqkPfT+D+BRBeA5L/gs8egrl94NR+b1cnIiIiIiLFnIK3SEFddxs8tgFuegZ8/GHvtzC9Nax9DbLSvV2diIiIiIgUUwreIp4wl4Obn4PH1kPNmyDrPKx+Bd5pC/tWe7s6EREREREphhS8RS5HpWuh3xLo/QEER8Hff8DHvWDhIDib4O3qRERERESkGFHwFrlcJhM0uhtGbIaWj4LJB3YtgqktYOO7YLN6u0IRERERESkGFLxFrlRAGHR7DR7+DqreAOnJ8PU/4L2OcGSrt6sTEREREREvU/AWuVqqNoMhq6D7G/YwnvAzvH8LfDUKzp32dnUiIiIiIuIlCt4iV5OPL7QYDCO2QON7AQO2zoa3m8OO/2rtbxERERGRMkjBW6QwBEfCXe/CgGVQqS6knYTFQ2FOd0jc7e3qRERERESkCCl4ixSmmBth6DqIewH8AuHgjzDjRlg5HjJSvV2diIiIiIgUAQVvkcLmZ4Ybn4QRm6BuN7BlwY9TYFor+G2Zt6sTEREREZFCpuAtUlTCr4H7/gv3/hfCoiHpMMy7H+beC6cPers6EREREREpJAreIkXt+m4wfKO9F9zHD37/2t77/cNkyMrwdnUiIiIiInKVKXiLeIM5yP7c99AfocaNkHUO4l+EGe1g//ferk5ERERERK4iBW8Rb4q8HgYshTvfhaAIOPk7fNgTPn8EUhK9XZ2IiIiIiFwFCt4i3mYyQZN7YcRmaD4YMMHP8+1rf2+aCTartysUEREREZEroOAtUlwEloceb8CQeKjSBNKTYPnT8H4cHN3u7epEREREROQyKXiLFDfVY+Hh1dD1dbCEwtFt8F4nWPY0nDvj7epERERERMRDCt4ixZGPL7R6BEZsgUb3AAZsnglTW8DPC8AwvF2hiIiIiIgUkIK3SHEWEgW934d+X0LFOpCaCJ8PgY9uhxO/e7s6EREREREpAAVvkZKgVkcY9j+4+TnwC7AvOfZOW/sSZBlp3q5OREREREQuQsFbpKTws8BNz8BjG+Da28CWCT9Mhumt4PcV3q5ORERERETyoeAtUtJUqAn3fwZ9P4HQanDmEMztA/MegDOHvV2diIiIiIjkoeAtUhKZTFCvJwzfBG0fBx8/+G0pTGsJP74F1kxvVygiIiIiIhcoeIuUZJZguO0lePQHuKYNZKbBynEwoz0c/J+3qxMRERERERS8RUqHqPowYDncMR3KVYQTu2F2V/hiGKSe9HZ1IiIiIiJlmoK3SGnh4wPNHrCv/X1Df/u+n+bC27GwZTbYbN6tT0RERESkjFLwFiltylWA2/8PBq+CqEZw/gwsHQUf3ArHfvJ2dSIiIiIiZY6Ct0hpFd0CHlkDnSeCORj+2gLvdYSvx8D5ZG9XJyIiIiJSZih4i5Rmvn7Q5jEYsRka3AmGDTa+A1NbwK5FYBjerlBEREREpNRT8BYpC0Krwj1z4MHPoUItSEmAhYPgk7vg733erk5EREREpFRT8BYpS+rcAsPWQ8ex4GuBfd/B9Naw+lXIPO/t6kRERERESqXLCt7Tpk0jJiaGgIAAWrVqxaZNmy7a/syZMwwfPpwqVapgsVi47rrrWL58+WUVLCJXyD8AOo6Bx9ZD7VvAmgFrJ9kD+B+rvF2diIiIiEip43Hwnj9/PqNHj2b8+PFs27aNJk2a0LlzZxITE922z8jI4NZbb+XAgQMsXLiQPXv2MHPmTKpVq3bFxYvIFahYGx5cZB+CHlIFTu+HT3rDZ/0g+ai3qxMRERERKTVMhuHZ7EqtWrWiRYsWTJ06FQCbzUZ0dDQjR45kzJgxLu1nzJjB66+/zm+//Ya/v/9lFZmcnExYWBhJSUmEhoZe1jlE5CLSz8LqibBxBhhW+yzoHcdCq6H2CdpERERERMSJJznVox7vjIwMtm7dSlxcXM4JfHyIi4tj/fr1bo9ZsmQJbdq0Yfjw4URFRdGwYUNeffVVrFarJ5cWkcJkCYEur8Kja6F6S8hIgW//Be91gEMbvV2diIiIiEiJ5lHwPnnyJFarlaioKKf9UVFRJCQkuD3mzz//ZOHChVitVpYvX87zzz/P5MmTefnll/O9Tnp6OsnJyU6biBSByo1g0Aro+X8QWB6O74JZt8GXIyDtlLerExEREREpkQp9VnObzUZkZCTvvfcesbGx9O3bl3/961/MmDEj32MmTpxIWFiYY4uOji7sMkUkm48PxPaHEVuh2YP2fds/hrdjYdvHYLN5tz4RERERkRLGo+BdqVIlfH19OX78uNP+48ePU7lyZbfHVKlSheuuuw5fX1/Hvnr16pGQkEBGRobbY8aOHUtSUpJjO3z4sCdlisjVEFQR7pgGA7+ByPpw7hQsGQGzu0LCLm9XJyIiIiJSYngUvM1mM7GxscTHxzv22Ww24uPjadOmjdtj2rVrxx9//IEtVy/Z77//TpUqVTCbzW6PsVgshIaGOm0i4iU12sCj38NtL4N/EBzeAO/eBCv+ZZ+UTURERERELsrjoeajR49m5syZfPjhh+zevZthw4aRmprKwIEDAejXrx9jx451tB82bBinTp3iiSee4Pfff2fZsmW8+uqrDB8+/Oq9CxEpXL7+0HYkjNgE9XraZz5fPxWmtoRfvwTPFkcQERERESlTPF4nqG/fvpw4cYJx48aRkJBA06ZN+eabbxwTrh06dAgfn5w8Hx0dzYoVK3jyySdp3Lgx1apV44knnuCf//zn1XsXIlI0wqpD30/g929h+dNw5qB93e86t0K316BCLW9XKCIiIiJS7Hi8jrc3aB1vkWIo8xz88Ab8OAWsGeAXAO2fgnZPgJ/F29WJiIiIiBSqQlvHW0TEwT8Qbv4XDPsf1OwAWedh9SswvQ3sW+3t6kREREREig0FbxG5MpWuhX5fQu8PIDgKTu2Dj3vBwkFwNsHb1YmIiIiIeJ2Ct4hcOZMJGt0NIzZDy0fB5AO7FsHbzWHDDLBmebtCERERERGvUfAWkasnIMw+ydrDq6FaLGSchW/+CTM7wZEt3q5ORERERMQrFLxF5Oqr2hQGr4Tub9jDeMLP8H4cfDUKzp32dnUiIiIiIkVKwVtECoePL7QYDCO2QpP7AAO2zrYPP98xV2t/i4iIiEiZoeAtIoUrOALunAEDlkGlupB2EhYPgzndIXG3t6sTERERESl0Ct4iUjRiboSh6yDuBfAvBwd/hBk3wspxkJHq7epERERERAqNgreIFB0/M9z4JAzfCHW7gy0LfnwLprWC35Z5uzoRERERkUKh4C0iRS/8GrhvLtw3D8KugaTDMO9+mHsvnD7o7epERERERK4qBW8R8Z66XWH4BrhxNPj4w+9f23u/v/8PZGV4uzoRERERkatCwVtEvMscBHHjYdiPENMess7Bdy/BjHaw/3tvVyciIiIicsUUvEWkeIioC/2/gjvfg6AIOPk7fNgTFj0MKYnerk5ERERE5LIpeItI8WEyQZO+MGIztBgCmGDnZ/a1vzfNBJvV2xWKiIiIiHhMwVtEip/A8tB9MjwcD1WaQnoSLH8a3r8F/trm7epERERERDyi4C0ixVe1WHj4O+j2H7CEwtHtMPNmWPY0nDvj7epERERERApEwVtEijcfX2j5MIzYAo3uAQzYPBOmtoCfPwPD8HaFIiIiIiIXpeAtIiVDSBT0fh/6LYGK10JqInz+sH0CthO/e7s6EREREZF8KXiLSMlSq4N96bGbnwO/ADjwA7zTFuJfhIw0b1cnIiIiIuJCwVtESh4/C9z0DAzfCNd2Blsm/DAZpreCPd94uzoREREREScK3iJScpWPgfvnQ99PIbQ6nDkE/+0L8x6AM4e9XZ2IiIiICKDgLSIlnckE9XrYe7/bPg4+fvDbUpjWEtZNAWumtysUERERkTJOwVtESgdLMNz2Ejz6A1zTFjLTYNV4mNEeDvzo7epEREREpAxT8BaR0iWqPgxcDndMh3IV4cRumNMNvhgGqSe9XZ2IiIiIlEEK3iJS+phM0OwB+9rfsQPs+36aC2/HwpZZYLN5tTwRERERKVsUvEWk9CpXAXq+BYNXQVQjOH8Glj4JH9wKx37ydnUiIiIiUkYoeItI6RfdAh5ZA13+DeYQ+GsLvNcRvv4nnE/2dnUiIiIiUsopeItI2eDrB62HwYjN0OAuMGywcQZMbQG7FoFheLtCERERESmlFLxFpGwJrQL3zIYHP4cKtSAlARYOgo/vhJN/eLs6ERERESmFFLxFpGyqcwsMWw8dnwVfC/y5Gt5pA9+9ApnnvF2diIiIiJQiCt4iUnb5B0DHf8Jj66H2LWDNgO9fg+ltYO8qb1cnIiIiIqWEgreISMXa8OAiuOdDCKkKp/fDp73hs36Q9Je3qxMRERGREk7BW0QE7Gt/N+gFIzZBmxFg8oVfv4RpLeF/U8Ga5e0KRURERKSEUvAWEcnNEgKdX4FH10L1lpCRAt/+C97rAIc2eLs6ERERESmBFLxFRNyp3AgGrYDb34bA8nB8F8zqDF+OgNS/vV2diIiIiJQgCt4iIvnx8YEb+sGIrdDsQfu+7R/D1Oaw7SOw2bxbn4iIiIiUCAreIiKXElQR7phm7wGPbADnTsGSkTC7CyTs8nZ1IiIiIlLMKXiLiBTUNa3tz37f9jL4B8HhjfDuTbDiX5B+1tvViYiIiEgxpeAtIuIJX39oOxJGbIZ6t4NhhfVTYWpL+GUxGIa3KxQRERGRYkbBW0TkcoRVg74fwwMLoXwMnD0KC/rDp3fDqT+9XZ2IiIiIFCMK3iIiV+LaW+GxDXDTP8DXDH+sgmmtYc0kyEr3dnUiIiIiUgwoeIuIXCn/QLj5XzBsPdTsANZ0WPMqTG8D+1Z7uzoRERER8TIFbxGRq6VSHej3JfT+AIKj4NQ++LgXLBgIyce8XZ2IiIiIeImCt4jI1WQyQaO77ZOvtRoKJh/45XOY2gI2vAPWLG9XKCIiIiJFTMFbRKQwBIRB10nw8GqoFgsZZ+GbMTCzExzZ4u3qRERERKQIXVbwnjZtGjExMQQEBNCqVSs2bdpUoOPmzZuHyWSiV69el3NZEZGSp2pTGLwSerxpD+MJP8P7cfDVE5B2ytvViYiIiEgR8Dh4z58/n9GjRzN+/Hi2bdtGkyZN6Ny5M4mJiRc97sCBAzz99NO0b9/+sosVESmRfHyh+SAYsRWa3A8YsHWOffj5jrla+1tERESklPM4eL/xxhs8/PDDDBw4kPr16zNjxgzKlSvHrFmz8j3GarXywAMPMGHCBGrVqnVFBYuIlFjBEXDnOzBgOURcD2knYfEwmN0NEnd7uzoRERERKSQeBe+MjAy2bt1KXFxczgl8fIiLi2P9+vX5Hvfiiy8SGRnJ4MGDC3Sd9PR0kpOTnTYRkVIjph08+gPETQD/cnDofzDjRlg5DjJSvV2diIiIiFxlHgXvkydPYrVaiYqKctofFRVFQkKC22PWrVvHBx98wMyZMwt8nYkTJxIWFubYoqOjPSlTRKT48zPDjaNg+Eao2x1sWfDjWzC1JexequHnIiIiIqVIoc5qfvbsWR566CFmzpxJpUqVCnzc2LFjSUpKcmyHDx8uxCpFRLwo/Bq4by7cN8/+7+QjMP8B+O+9cPqAt6sTERERkavAz5PGlSpVwtfXl+PHjzvtP378OJUrV3Zpv2/fPg4cOEDPnj0d+2w2m/3Cfn7s2bOH2rVruxxnsViwWCyelCYiUrLV7Qo1O8AP/4Ef/w9+/wb+XAs3PQ1tH7f3kIuIiIhIieRRj7fZbCY2Npb4+HjHPpvNRnx8PG3atHFpf/3117Nz50527Njh2G6//XY6derEjh07NIRcRCQ3czm4ZRwM+xFi2kPWOfjuJZjRDvZ/7+3qREREROQyedTjDTB69Gj69+9P8+bNadmyJVOmTCE1NZWBAwcC0K9fP6pVq8bEiRMJCAigYcOGTseHh4cDuOwXEZELIupC/69g5wJY8Syc/B0+7AmN+sBtL0NI1KXPISIiIiLFhsfBu2/fvpw4cYJx48aRkJBA06ZN+eabbxwTrh06dAgfn0J9dFxEpPQzmaBxH7j2Nnuv9+YPYOdn8PsKuOV5+7rgPr7erlJERERECsBkGMV/6tzk5GTCwsJISkoiNDTU2+WIiBS9v7bC0tFwbIf98ypNocebUO0Gb1YlIiIiUmZ5klMVvEVESgqbFbbMgviXID0JMEHMjRDdEqq3hOotIKiit6sUERERKRMUvEVESrOzx+Hb5+xDz/OqWMcewqNb2D9G1tOQdBEREZFCoOAtIlIWnPgdDv4IRzbD4U3w917XNuYQ+3B0R694cyhXoehrFRERESllFLxFRMqitFNwZAsc2WQP4n9thYwU13YVr70QxFvYP0Zcr15xEREREQ8peIuIiP2Z8MTdF4L4ZvvHv/9wbWcJhWqxuXrFYyGwfNHXKyIiIlKCKHiLiIh7qX/bh6Y7esW3QWaqa7tKdXOeE49uaf9cS0WKiIiIOCh4i4hIwVizIPFX517xU3+6trOE2XvCsyduq9YcAsOLvFwRERGR4kLBW0RELl/qyZwJ245stj8rnpmWp5EJIurmPCdevSVUuk694iIiIlJmKHiLiMjVY82CxF/sQfzwJnuv+OkDru0CwuxBPHeveIB+ZouIiEjppOAtIiKFK+VEznPiRzbbnxXPOpenkcm+jnjuXvGKddQrLiIiIqWCgreIiBQtayYc35XznPjhTXDmoGu7gPBcQbyFfTZ19YqLiIhICaTgLSIi3nf2eK4Z1DfD0W2QdT5PIxNE1rcHcUeveG0wmbxSsoiIiEhBKXiLiEjxk5UBx3fm6hXfDEmHXNsFVrjQK37hefFqsWAJLvp6RURERC5CwVtEREqGswk5E7Yd3gxHt4M13bmNyQciGzivK16hlnrFRURExKsUvEVEpGTKyoCEnc4TtyUddm1XruKFGdRbQHQrqHYDmIOKvl4REREpsxS8RUSk9Eg+mhPCD2+CYzvAmuHcxuQLUQ1ynhOPbgHla6pXXERERAqNgreIiJReWelw7GfnXvHkv1zbBUXk6hVvCVVvAHO5oq9XRERESiUFbxERKVuS/sp5TvzIJjj2k/te8coNc54Tr94CyseoV1xEREQui4K3iIiUbZnnIeHnXBO3bYKzx1zbBUXmhPDollC1GfgHFn29IiIiUuIoeIuIiORmGJB0JE+v+M9gy3Ru5+MHlRs594qHX6NecREREXGh4C0iInIpmeftE7XlXs4sJcG1XXBUTo949ZZQtal6xUVERETBW0RExGOGYV+6LPcM6gk/gy3LuZ2Pv71XPLplThgPq65ecRERkTJGwVtERORqyDwHR7c7h/HURNd2IVWce8WrNAH/gKKvV0RERIqMgreIiEhhMAw4czDnOfHDmyBhJxhW53Y+/vbwnXvitrDq3qlZRERECoWCt4iISFHJSLP3iueeuC31hGu7kKoQ3eLCxG2toEpj8LMUfb0iIiJyVSh4i4iIeIthwOkDOUPTD2+E47+49or7mqFKU+de8dCq3qhYRERELoOCt4iISHGSkQp/bXPuFU/727VdaPVcveItoXJj8DMXfb0iIiJySQreIiIixZlhwKk/c3rFj2y60Ctuc27na7EvX5Z74rbQKl4pWURERJwpeIuIiJQ06SlwdJvzDOrnTrm2C4u+EMRb2XvHoxqpV1xERMQLFLxFRERKuuxe8cMbc8J44q+uveJ+AVC1mXOveEiUd2oWEREpQxS8RURESqP0s/DX1pznxI9shnOnXduFX5PznHj1FlC5Efj6F329IiIipZiCt4iISFlgGPD3HznPiR++0CtOnv+1+wXae8VzT9wWHOmVkkVEREoLBW8REZGy6nzyhV7xTTm94ueTXNuVj3HuFY9qCL5+RV6uiIhISaXgLSIiInY2G/y917lX/MRvuPSK+5eDqjc494oHVfJKySIiIiWBgreIiIjk73wSHNmSazmzLZDurle8Zk6PeHRLiGygXnEREZELFLxFRESk4Gw2OPn7hR7xC8PTT/zm2s4/CKrdkDN7evUWEFSx6OsVEREpBhS8RURE5MqcOw1HtuaE8b+2Qnqya7sKtfP0itcHH9+ir1dERKSIKXiLiIjI1WWzwok9Oc+JH9lk7yXPyxxs7xXPPXFbuQpFX6+IiEghU/AWERGRwpd2Ks8M6lsh46xru4p1LgTxCxO3RdZTr7iIiJR4Ct4iIiJS9GxW+7Phhzfm9Ir//YdrO3MIVI/N1SveHALLF329IiIiV0DBW0RERIqHtFO5Zk/fBH9tg4wU13aVrnPuFY+4Hnx8ir5eERGRAlLwFhERkeLJZoXEX3NmTz+8CU7tc21nCYVqsblmUG8OgeFFXq6IiEh+FLxFRESk5Ej92x7Cc8+gnpnm2q5SXXsQzw7jla5Tr7iIiHiNgreIiIiUXNYsSPzFuVf89H7XdgFhUK15zuzp1Zvb94mIiBSBQg/e06ZN4/XXXychIYEmTZrw9ttv07JlS7dtZ86cyUcffcSuXbsAiI2N5dVXX823vTsK3iIiImVcyolcveKb4eg2N73iJvuz4dnPiUe3hIrXqldcREQKRaEG7/nz59OvXz9mzJhBq1atmDJlCgsWLGDPnj1ERka6tH/ggQdo164dbdu2JSAggEmTJvHFF1/wyy+/UK1atav+hkRERKQMsGbB8V3OE7edPuDaLiDc3hOePXFbteYQoN8lRETkyhVq8G7VqhUtWrRg6tSpANhsNqKjoxk5ciRjxoy55PFWq5Xy5cszdepU+vXrV6BrKniLiIjIJaUk5oTww5vh6HbIOpenkQki6+fpFa8DJpNXShYRkZLLk5zq58mJMzIy2Lp1K2PHjnXs8/HxIS4ujvXr1xfoHGlpaWRmZlKhQgVPLi0iIiJyccGRUK+HfQOwZkLCTude8TOH7M+PJ/4CW+fY2wWWv/CMeHaveCxYQrz2NkREpPTxKHifPHkSq9VKVFSU0/6oqCh+++23Ap3jn//8J1WrViUuLi7fNunp6aSnpzs+T05O9qRMEREREfD1h2o32LdWj9r3nT2eM3v6kQu94udOw95v7RuAycfeK169Rc4M6hVrq1dcREQum0fB+0r9+9//Zt68eaxZs4aAgIB8202cOJEJEyYUYWUiIiJSJoREQb2e9g0gKwOO77QPTT+80R7Gkw7bnx8/vgu2zra3K1fxQq/4hTBe9QawBHvvfYiISIni0TPeGRkZlCtXjoULF9KrVy/H/v79+3PmzBm+/PLLfI/9z3/+w8svv8yqVato3rz5Ra/jrsc7Ojpaz3iLiIhI4Us+lqdXfAdY053bmHwgqkHOc+LVW0CFWuoVFxEpQwrtGW+z2UxsbCzx8fGO4G2z2YiPj2fEiBH5Hvfaa6/xyiuvsGLFikuGbgCLxYLFYvGkNBEREZGrI7QK1L/DvgFkpdufFc89cVvyEfu+hJ2w5QN7u3KVLvSIX3hevNoNYA7y3vsQEZFiw+Oh5qNHj6Z///40b96cli1bMmXKFFJTUxk4cCAA/fr1o1q1akycOBGASZMmMW7cOObOnUtMTAwJCQkABAcHExysIVoiIiJSzPlZLixJ1hx4zL4v+WhOj/jhTXBsB6SdhN+/tm8AJl97r3h0S4huZQ/l5WPUKy4iUgZ5HLz79u3LiRMnGDduHAkJCTRt2pRvvvnGMeHaoUOH8PHxcbR/5513yMjI4O6773Y6z/jx43nhhReurHoRERERbwitCg162Tew94of+8m5V/zsUUj42b5tft/eLigiZ/b06i2hajMwl/PWuxARkSLi8Tre3qB1vEVERKTESTqSp1f8J7BlOrfx8YOohjmzp0e3gPAa6hUXESkBPMmpCt4iIiIiRSHzvD1855647ewx13ZBkTkTtkVf6BX3Dyz6ekVE5KIUvEVERESKO8OwL12Wu1c84WewZTm38/GDyo1ynhOPbglh0eoVFxHxMgVvERERkZIo85x9+bLcveIpx13bBVfOeU48uiVUaQr+AUVdrYhImabgLSIiIlIaGAacOZTTI35kk30JM5decX+o0th54raw6uoVFxEpRAreIiIiIqVVRpp9+TLHEPWNkHrCtV1IlQtD01td6BVvYl8aTURErgoFbxEREZGywjDg9IE8veK7wLA6t/M128O3U694Na+ULCJSGih4i4iIiJRlGalwdLvzxG1pJ13bhVbLmbCtekv7cHX1iouIFIiCt4iIiIjkMAw4vR8Ob86ZuO34L256xS32XvHcy5mFVvVOzSKSL8MwsNoMrNkfc2/u9l3Yn2U1sBkGWTYDW97XLuxz+njhGKvhui/7PO6uk93WqSarZ/UO7ViLm6+P8vaX+qI8yal+RVSTiIiIiHiLyQQVatm3Jn3t+9JT7L3i2UH88CY4d8r++ZFNOceGVrcPTY+sD4Hlc23hOf+2hIGPj1fempQuFw2U+YVHN/tyf3QX+nIHyvyuU5B9eUOqy2tONdmw2SDLZsNqkKdOGzYD5zY2sNpsF+p1fs1W7LtOr9ydyaXrURj1eIuIiIiIvVf81J85z4kf3gyJv4BhK8DBJnsQDwjPE87dhPTsLbutn7lQ31ZxkR0oXcOYc+9g3tcuGQQvER7d9WzmFxqzw6vb13IFREcgdIRH50BpzRMac3pBXQNl3prKQqAsbCYT+JpM+Po4b34+JnxMFz7mfs2Up03utr72j75u9uVum/tcPnnOYz+/D74+OD5mH28/zvk1Xx8ffE0mGlUL45qK5bz95bwo9XiLiIiIiGdMJqhY2741vc++L/0s/LXNHsTPHIJzp+HcmQsfL/w7MxUwcvad3u/RZa1+5ciyhJNlDiPLHEqGOZxM/zAy/MNI9w8lwxzGeb9Q0v1COecXSrpfGOf8QjlvCrjQC+h+WGx2+HMeFps79OWEx5wQ6G6f+97Wi4Vfd72gCpRXzscEvj4XD32+vs5BMjsougui9oBHPvucw6Cfj8+FAOoaEN0FUZ8L9eTdlzv8+uYTcv1c6synfjfvNzv4SvGj4C0iIiJSQhmGQabV4HyWlfOZVtIzbZzPtHI+0+bYd96xz8r5LBvpmXn2Z+VuYyM973FZVtIzm5FpbeK2F9TfyCSMVMJMKYSTQpgplXBSCTelXNiXemFfzufhphRCScPHZOCblYZvVhqW1KMevfdMw5czBJFkBHOGYM4YQSRlf24Ec5YgzhjBJF/4eCbX5zaK17B4H9OFYOeTt6cyV+hzhL/sz3OFQFOeNrnCnK+bfW5D66UCbUGCYH5hsAA9qy773IRUk9allxJMwVtERETkKrHajAvBNVfYzRWC0/MJuzkB1/1x2W3Ss3Jez/63t3tSM/DnBOGcMMLtO3LV4whbpuxgB36+9p5Df5NBmM85wk0phJsuBHUjhVBSCSWFUM4SaqQQYpwl2Egh2JZCsO0sQbaz+JOJv8lKBMlEmJI9rvm8XwgZfqH2HnX/MDLM9h72LHMYmeYwsszhWAPsH22WMKwB4RgB5cE/INfQWfc9lh6HVAVKkTJBwVtERERKJcMwyLDa7L24l+gFtofi/HuB7aH4Ir3EF/ZlWr2Xgk0mCPDzJcDfhwB/XwL8fbH4Zf/7wscLr1tytbNkv+7n69z2wr7s182+Phd6JXP3rHJhOK3zvuye2UIJlIYBmWl5hrzn2s67259k/5hxFoCArLMEZJ2F8395dm2/QDfPrYe7Preee7OEgyXU/g0SkTJLwVtERESKRJbVljPUOcu5Zzf9Ir3AOcE5n/Ccd/h0Vk6PsDenkDX7+mDJE2JzB1uLmxCct53lQoAO8MsJ0y4hOlcwLhM9pyYTmIPsW5iHsx5bM3MCu9uA7mZfdjvDBlnn4Ow5OOvZsHhMvpcO6O4moQsIB1/9ui5SGui/ZBERkTLIMIxcw5ZtrsOjcw9pdttTfKnh067nzfLimGgfEznB1c9dT2/eAOxrD8359AI7Xs8VgHMHY4ufL76a4Kj48fWH4Aj75gmbzd5bftGQfsZ9mM86b18vPe1v++Ypc4ibUB7uPrjnDvT+geplFylGFLxFRESKgUyrc09vdihOv0gvsGN4dJ7eY7fDp908K+xNZj+ffHtxLZfoBXaE5ov0AucNxv6+eo5WroCPDwSE2bfyMZ4dm3nONaS77W3PHebPQHqS/fiMs/Yt6ZBn1/W1XGIpt7z7w7Umu0ghUvAWERHJw2bL1RuclU+PsLvZn/MEXNfh07mfF3buXbZ6sTfY18fkFGItjmd7ffLtzXUOu3lez6cX2NFT7Oej5W6k7PAPtG+hVTw7zpoF55MKPiQ+d6C3ZYE1HVIS7JtHTAUYFu9mX0B4mVmTXeRyKHiLiEixdrHlkpzCcJ6Am36RXuD8h0/b22R4uTc474RYjs/9cvXiOgVed5Np5T9hVt7g7O+r3i2RYsfXD4Iq2jdPGAZkpBTsufVLrcnuKf+gPOE8/OLD4bM3c5CGxUupp+AtIiIesdmMfHtxswNvupsw6zzMOb/jc02SlSske3O5JD8fU57A6voscN7Zn12f+3WeTOtiwdjsq95gEbkCJhNYQuxb+DWeHZuVXsBh8XnbJAGGPbhnpkLyEc+u6+Nf8JDu1BsfBj6+nl1LxEsUvEVESrCCLJfk1PObd/Znd8siucwYnft5YhsZVu/2Bud9Fth5Zmf3vcABfrleL2AvcPZxfuoNFpGyws8CIVH2zRM2m/2Z9PyeV79YkLdmgC0TUk/YN08FhHk+JD7Qvia7SFFS8BYRKWayrDZOpWaQeDadE9lbSq5/5/o8NSPLq8sl+fua8pkMqyCzP19q+LS7CbfKyHJJIiIliY9PTrD1hGNN9ksNic9/TXb7c/BJntd8qTXZ812XPUTD4uWyKHiLiBQBwzBIPp/lEpwTz553CtQnU9L5OzXjssK0yYTTs8DuZ3/O5/WLLJtkcdMLnP0ssXqDRUTksjmtyV7ds2MvuSb7RcL81VyTvaDrsmtN9jJP330RkSuQnmV12xOd+/PEZPtHTybs8jFBxWALEcEWIkLsW2RIzr8jgi1UCrEQYvFzBGOzr3qDRUSkjCiJa7JbQi++lJvWZC/VFLxFRPKw2QxOp2Xk9Eon5z/UO+lcpkfnDgnwcwTniDxBOjI0wLG/QpAZX02wJSIicnVd7TXZLzkB3ZmcNdnTk+3bZa/J7i6kh+f/fLvWZC9WFLxFpMxITc9y3yudPeQ7JXu4d4ZHayr7+5pyBekAl0Cd3VtdKdhCoFmzr4qIiJRIJW1NdpNPzuRzBR0SrzXZC42Ct4iUaFlWG3+nZrh9Xjp3wE48m05ahtWjc1cIMrv0TEfmCdQRIRbCAv01xFtERETc89aa7IbtKq7Jns+Q+LyBXmuy50vBW0SKHcMwSD6XxYmU884ze7sJ1KfSPJuILNDfl8jQPEO9XcJ1ABWDzfhr4jARERHxlhK7JnsBnlvPG+rLwJrsCt4iUmTOZ1ovujRW7s2TtaJ9TFDJ3SRkwa5Dv4PMvuqdFhERkdLNq2uyJ9o3T+Vdk73VUKjbxfPzFFMK3iJyRWw2g1NpGU5Dup0Ddc7w7+TzWR6dOzR7IrLsZ6eD8wTrC1v5cpqITEREROSKFac12RvdfbXeVbGg4C0ibqWkZ+XphT6fa+3pnP1/p3o2EZnZ18d58rF8ZviOCLEQ4F+6hxyJiIiIlAqFsSZ7dMvCqNRrFLxFypBMq42/U9xMRJbiOuzb04nIKgaZnYNznueoI0MsRAQHEBrop6HeIiIiImJ3uWuylzAK3iIlnGEYJJ3LdAnOLpOSpaRzKjXDo3MHmX3z7ZmOzPXsdIUgTUQmIiIiIpIfBW+RYip7IrJEd73SZ3OenT6ZkuHRRGS+PiYqBZtzwnPeYd65AnaQRT8iRERERESulH6rFilCVpvBqQtrTp9ISScx+Xy+M3yf9XAisrBAfze90q6Bunw5Mz6aiExEREREpMgoeItcIcMwnCciyxOkc/dY/52SjgfzkGH288mzNJZrr3RkaACVgs1Y/DQRmYiIiIhIcaTgLZKPjCwbf6deCM/JF197+lxmwSciM5nsE5FVCnZ9VjpvwA4N0ERkIiIiIiIlnYK3lCmGYXAmLdNNr/R5lzB9Oi3To3MHW/zc90znCtSRFyYi89NEZCIiIiIiZYaCt5QK5zKsF4Lzebe90tnDvU+mpJNpLfhYbz8fE5WCLUSG5j/UO3srZ9Z/TiIiIiIi4kpJQYotq81wDPV2WR7rQqA+eeHzs+meTUQWXs4//0nIgnOGfocH+msiMhERERERuSIK3lKkDMPgbO6JyC6y9vSpVM8mIrP4+bj2TOcK0dnhuqImIhMRERERkSKk4C1XRUaWjZMpeXql8w79vhCuz2cWfM1p+0Rk7nqlXZ+hDrFoIjIRERERESl+FLwlXzabwZlzmS4h2mWG75R0zng4EVnIhYnIKuWeeMzNc9QVymkiMhERERERKdkUvMugtIwst8O8c3+emGyfiCzLg7He/r6mS05AFhkSQKVgC4FmDfUWEREREZGy4bKC97Rp03j99ddJSEigSZMmvP3227Rs2TLf9gsWLOD555/nwIEDXHvttUyaNIlu3bpddtHiKstq41Rqhn2od0o6Jy6y7nSKhxORlS/n7xSc3Q71DrYQponIREREREREXHgcvOfPn8/o0aOZMWMGrVq1YsqUKXTu3Jk9e/YQGRnp0v5///sf9913HxMnTqRHjx7MnTuXXr16sW3bNho2bHhV3kRpZRgGyeez3PZMO9advrBE1t+pGRgeTEQW4O+Tb4jO/Sx1xSALZj8N9RYREREREblcJsPwJK5Bq1ataNGiBVOnTgXAZrMRHR3NyJEjGTNmjEv7vn37kpqaytKlSx37WrduTdOmTZkxY0aBrpmcnExYWBhJSUmEhoZ6Um6xlJ5l5WRKBonJ5/Md6p397/Ssgk9E5mOCisG5gnM+w74jQwMIMvtqIjIREREREZHL5ElO9ajHOyMjg61btzJ27FjHPh8fH+Li4li/fr3bY9avX8/o0aOd9nXu3JnFixd7cukS4UxaBgnJ551CtLu1p5POeTgRWYCfS5B29Fbn2l8hyIyvhnqLiIiIiIgUKx4F75MnT2K1WomKinLaHxUVxW+//eb2mISEBLftExIS8r1Oeno66enpjs+Tk5M9KdNrRs3fwZo9JwrU1uzrkzOrt5ue6dwzfAf4ayIyERERERGRkqpYzmo+ceJEJkyY4O0yPBZ5odfZMaQ7xP1Q74gQ+0RkGuotIiIiIiJS+nkUvCtVqoSvry/Hjx932n/8+HEqV67s9pjKlSt71B5g7NixTsPTk5OTiY6O9qRUr5jUu7HCtIiIiIiIiDjxaLpqs9lMbGws8fHxjn02m434+HjatGnj9pg2bdo4tQdYuXJlvu0BLBYLoaGhTltJoNAtIiIiIiIieXk81Hz06NH079+f5s2b07JlS6ZMmUJqaioDBw4EoF+/flSrVo2JEycC8MQTT9ChQwcmT55M9+7dmTdvHlu2bOG99967uu9EREREREREpBjyOHj37duXEydOMG7cOBISEmjatCnffPONYwK1Q4cO4eOT05Hetm1b5s6dy3PPPcezzz7Ltddey+LFi7WGt4iIiIiIiJQJHq/j7Q2lbR1vERERERERKdk8yakePeMtIiIiIiIiIp5R8BYREREREREpRAreIiIiIiIiIoVIwVtERERERESkECl4i4iIiIiIiBQiBW8RERERERGRQqTgLSIiIiIiIlKI/LxdQEFkLzWenJzs5UpEREREREREcvJpdl69mBIRvM+ePQtAdHS0lysRERERERERyXH27FnCwsIu2sZkFCSee5nNZuPo0aOEhIRgMpm8XU6+kpOTiY6O5vDhw4SGhnq7HClBdO/I5dB9I5dD941cDt03cjl038jlKEn3jWEYnD17lqpVq+Ljc/GnuEtEj7ePjw/Vq1f3dhkFFhoaWuxvEimedO/I5dB9I5dD941cDt03cjl038jlKCn3zaV6urNpcjURERERERGRQqTgLSIiIiIiIlKIFLyvIovFwvjx47FYLN4uRUoY3TtyOXTfyOXQfSOXQ/eNXA7dN3I5Sut9UyImVxMREREREREpqdTjLSIiIiIiIlKIFLxFRERERERECpGCt4iIiIiIiEghUvD20LRp04iJiSEgIIBWrVqxadOmi7ZfsGAB119/PQEBATRq1Ijly5cXUaVSnHhy38yZMweTyeS0BQQEFGG1Uhx8//339OzZk6pVq2IymVi8ePElj1mzZg033HADFouFOnXqMGfOnEKvU4oXT++bNWvWuPy8MZlMJCQkFE3BUixMnDiRFi1aEBISQmRkJL169WLPnj2XPE6/45Rtl3Pf6Hcceeedd2jcuLFjje42bdrw9ddfX/SY0vKzRsHbA/Pnz2f06NGMHz+ebdu20aRJEzp37kxiYqLb9v/73/+47777GDx4MNu3b6dXr1706tWLXbt2FXHl4k2e3jcAoaGhHDt2zLEdPHiwCCuW4iA1NZUmTZowbdq0ArXfv38/3bt3p1OnTuzYsYNRo0YxZMgQVqxYUciVSnHi6X2Tbc+ePU4/cyIjIwupQimO1q5dy/Dhw9mwYQMrV64kMzOT2267jdTU1HyP0e84cjn3Deh3nLKuevXq/Pvf/2br1q1s2bKFm2++mTvuuINffvnFbftS9bPGkAJr2bKlMXz4cMfnVqvVqFq1qjFx4kS37fv06WN0797daV+rVq2MRx99tFDrlOLF0/tm9uzZRlhYWBFVJyUBYHzxxRcXbfOPf/zDaNCggdO+vn37Gp07dy7EyqQ4K8h9s3r1agMwTp8+XSQ1ScmQmJhoAMbatWvzbaPfcSSvgtw3+h1H3Clfvrzx/vvvu32tNP2sUY93AWVkZLB161bi4uIc+3x8fIiLi2P9+vVuj1m/fr1Te4DOnTvn215Kn8u5bwBSUlKoUaMG0dHRF/0roEg2/byRK9G0aVOqVKnCrbfeyo8//ujtcsTLkpKSAKhQoUK+bfQzR/IqyH0D+h1HclitVubNm0dqaipt2rRx26Y0/axR8C6gkydPYrVaiYqKctofFRWV77NwCQkJHrWX0udy7pu6desya9YsvvzySz755BNsNhtt27blyJEjRVGylFD5/bxJTk7m3LlzXqpKirsqVaowY8YMFi1axKJFi4iOjqZjx45s27bN26WJl9hsNkaNGkW7du1o2LBhvu30O47kVtD7Rr/jCMDOnTsJDg7GYrEwdOhQvvjiC+rXr++2bWn6WePn7QJExFmbNm2c/urXtm1b6tWrx7vvvstLL73kxcpEpLSpW7cudevWdXzetm1b9u3bx5tvvsnHH3/sxcrEW4YPH86uXbtYt26dt0uREqSg941+xxGw/79nx44dJCUlsXDhQvr378/atWvzDd+lhXq8C6hSpUr4+vpy/Phxp/3Hjx+ncuXKbo+pXLmyR+2l9Lmc+yYvf39/mjVrxh9//FEYJUopkd/Pm9DQUAIDA71UlZRELVu21M+bMmrEiBEsXbqU1atXU7169Yu21e84ks2T+yYv/Y5TNpnNZurUqUNsbCwTJ06kSZMmvPXWW27blqafNQreBWQ2m4mNjSU+Pt6xz2azER8fn+8zCW3atHFqD7By5cp820vpczn3TV5Wq5WdO3dSpUqVwipTSgH9vJGrZceOHfp5U8YYhsGIESP44osv+O6776hZs+Ylj9HPHLmc+yYv/Y4jYP/dOD093e1rpepnjbdndytJ5s2bZ1gsFmPOnDnGr7/+ajzyyCNGeHi4kZCQYBiGYTz00EPGmDFjHO1//PFHw8/Pz/jPf/5j7N692xg/frzh7+9v7Ny501tvQbzA0/tmwoQJxooVK4x9+/YZW7duNe69914jICDA+OWXX7z1FsQLzp49a2zfvt3Yvn27ARhvvPGGsX37duPgwYOGYRjGmDFjjIceesjR/s8//zTKlStnPPPMM8bu3buNadOmGb6+vsY333zjrbcgXuDpffPmm28aixcvNvbu3Wvs3LnTeOKJJwwfHx9j1apV3noL4gXDhg0zwsLCjDVr1hjHjh1zbGlpaY42+h1H8rqc+0a/48iYMWOMtWvXGvv37zd+/vlnY8yYMYbJZDK+/fZbwzBK988aBW8Pvf3228Y111xjmM1mo2XLlsaGDRscr3Xo0MHo37+/U/vPPvvMuO666wyz2Ww0aNDAWLZsWRFXLMWBJ/fNqFGjHG2joqKMbt26Gdu2bfNC1eJN2cs85d2y75X+/fsbHTp0cDmmadOmhtlsNmrVqmXMnj27yOsW7/L0vpk0aZJRu3ZtIyAgwKhQoYLRsWNH47vvvvNO8eI17u4ZwOlniH7Hkbwu577R7zgyaNAgo0aNGobZbDYiIiKMW265xRG6DaN0/6wxGYZhFF3/uoiIiIiIiEjZome8RURERERERAqRgreIiIiIiIhIIVLwFhERERERESlECt4iIiIiIiIihUjBW0RERERERKQQKXiLiIiIiIiIFCIFbxEREREREZFCpOAtIiIiIiIiUogUvEVERIqRAQMG0KtXL2+XISIiIleRn7cLEBERKStMJtNFXx8/fjxvvfUWhmEUUUUFs2bNGjp16sTp06cJDw/3djkiIiIljoK3iIhIETl27Jjj3/Pnz2fcuHHs2bPHsS84OJjg4GBvlCYiIiKFSEPNRUREikjlypUdW1hYGCaTyWlfcHCwy1Dzjh07MnLkSEaNGkX58uWJiopi5syZpKamMnDgQEJCQqhTpw5ff/2107V27dpF165dCQ4OJioqioceeoiTJ0/mW9vBgwfp2bMn5cuXJygoiAYNGrB8+XIOHDhAp06dAChfvjwmk4kBAwYAYLPZmDhxIjVr1iQwMJAmTZqwcOFCxznXrFmDyWRi2bJlNG7cmICAAFq3bs2uXbuu3hdVRESkBFDwFhERKeY+/PBDKlWqxKZNmxg5ciTDhg3jnnvuoW3btmzbto3bbruNhx56iLS0NADOnDnDzTffTLNmzdiyZQvffPMNx48fp0+fPvleY/jw4aSnp/P999+zc+dOJk2aRHBwMNHR0SxatAiAPXv2cOzYMd566y0AJk6cyEcffcSMGTP45ZdfePLJJ3nwwQdZu3at07mfeeYZJk+ezObNm4mIiKBnz55kZmYW0ldLRESk+DEZxe1BMhERkTJgzpw5jBo1ijNnzjjtHzBgAGfOnGHx4sWAvcfbarXyww8/AGC1WgkLC+Ouu+7io48+AiAhIYEqVaqwfv16Wrduzcsvv8wPP/zAihUrHOc9cuQI0dHR7Nmzh+uuu86lnsaNG9O7d2/Gjx/v8pq7Z7zT09OpUKECq1atok2bNo62Q4YMIS0tjblz5zqOmzdvHn379gXg1KlTVK9enTlz5lz0DwEiIiKliZ7xFhERKeYaN27s+Levry8VK1akUaNGjn1RUVEAJCYmAvDTTz+xevVqt8+L79u3z23wfvzxxxk2bBjffvstcXFx9O7d2+m6ef3xxx+kpaVx6623Ou3PyMigWbNmTvtyB/MKFSpQt25ddu/efbG3LCIiUqooeIuIiBRz/v7+Tp+bTCanfdmzpdtsNgBSUlLo2bMnkyZNcjlXlSpV3F5jyJAhdO7cmWXLlvHtt98yceJEJk+ezMiRI922T0lJAWDZsmVUq1bN6TWLxVLAdyYiIlI2KHiLiIiUMjfccAOLFi0iJiYGP7+C/68+OjqaoUOHMnToUMaOHcvMmTMZOXIkZrMZsA9zz1a/fn0sFguHDh2iQ4cOFz3vhg0buOaaawA4ffo0v//+O/Xq1buMdyYiIlIyaXI1ERGRUmb48OGcOnWK++67j82bN7Nv3z5WrFjBwIEDncJzbqNGjWLFihXs37+fbdu2sXr1akc4rlGjBiaTiaVLl3LixAlSUlIICQnh6aef5sknn+TDDz9k3759bNu2jbfffpsPP/zQ6dwvvvgi8fHx7Nq1iwEDBlCpUiWnmdtFRERKOwVvERGRUqZq1ar8+OOPWK1WbrvtNho1asSoUaMIDw/Hx8f9//qtVivDhw+nXr16dOnSheuuu47p06cDUK1aNSZMmMCYMWOIiopixIgRALz00ks8//zzTJw40XHcsmXLqFmzptO5//3vf/PEE08QGxtLQkICX331laMXXUREpCzQrOYiIiJSKNzNhi4iIlIWqcdbREREREREpBApeIuIiIiIiIgUIg01FxERERERESlE6vEWERERERERKUQK3iIiIiIiIiKFSMFbREREREREpBApeIuIiIiIiIgUIgVvERERERERkUKk4C0iIiIiIiJSiBS8RURERERERAqRgreIiIiIiIhIIVLwFhERERERESlE/w9GxMPNFGsSfQAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 1000x400 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "States:\n",
            " [[0.         1.        ]\n",
            " [0.09685631 0.37126203]\n",
            " [0.13283732 0.14222434]\n",
            " [0.14670236 0.07074541]]\n",
            "Controls:\n",
            " [[-0.62873797]\n",
            " [-0.22903769]\n",
            " [-0.07147893]]\n"
          ]
        }
      ],
      "source": [
        "# Solve the LQR problem with the classical Riccati equation method\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.linalg import solve_discrete_are\n",
        "\n",
        "# System definition\n",
        "A = np.array([[1.0, 0.1],\n",
        "        [0.0, 1.0]])\n",
        "B = np.array([[0.005],\n",
        "        [1.0]])\n",
        "Q = np.eye(2)\n",
        "R = np.array([[1.0]])\n",
        "x0 = np.array([0.0, 1.0])\n",
        "N = 3  # Input horizon length (state horizon would be N+1=4)\n",
        "\n",
        "# Backward Riccati recursion to compute K_t\n",
        "P = Q.copy()  # Terminal cost\n",
        "K_list = []\n",
        "\n",
        "for _ in range(N):\n",
        "    S = R + B.T @ P @ B\n",
        "    K = np.linalg.solve(S, B.T @ P @ A)  # K = (R + Bᵀ P B)⁻¹ Bᵀ P A\n",
        "    K_list.insert(0, K)  # prepend to use forward later\n",
        "    P = Q + A.T @ P @ A - A.T @ P @ B @ K  # update P for next step\n",
        "\n",
        "# Forward rollout of trajectory\n",
        "x = x0.reshape(-1, 1)\n",
        "states = [x.flatten()]\n",
        "controls = []\n",
        "\n",
        "for K in K_list:\n",
        "    u = -K @ x\n",
        "    controls.append(u.flatten())\n",
        "    x = A @ x + B @ u\n",
        "    states.append(x.flatten())\n",
        "\n",
        "states = np.array(states)\n",
        "controls = np.array(controls)\n",
        "\n",
        "# Plot\n",
        "plt.figure(figsize=(10, 4))\n",
        "plt.plot(states[:, 0], label='Position')\n",
        "plt.plot(states[:, 1], label='Velocity')\n",
        "plt.title('Riccat-equation-solved State Trajectories')\n",
        "plt.xlabel('Time step')\n",
        "plt.legend()\n",
        "\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# Optional: print values\n",
        "print(\"States:\\n\", states)\n",
        "print(\"Controls:\\n\", controls)\n"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "gpuType": "A100",
      "machine_shape": "hm",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
